“Calhoun 


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1990-06 


An analysis of MLAYER: a multilayer 
tropospheric propagation program 


Yeoh, Lean-Weng 


Monterey, California. Naval Postgraduate School 


NPS-62-90-009 
http://hdl.handle.net/10945/44295 


Copyright is reserved by copyright owner 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 


\§ D U DL EY research materials and institutional publications created by the NPS community. 
«iis Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed -- and published -- scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 


http://www.nps.edu/library Monterey, California USA 93943 


“aC alhoun 


Institutional Archive of the Naval Postgraduate School 








Author(s) Y eoh, Lean-Weng 


Anandysis of MLAY ER: amultilayer tropospheric propagation program 





Publisher Monterey, California. Naval Postgraduate School 





Issue Date 1990-06 





URL http://hdl handle. net/10945/44295 








This document was downloaded on J anuary 15, 2015 at 10:37:41 


th D U DLEY Calhoun is a project of the Dudley Knox Library at NPS, furthering the precepts and 
«iis goals of open government and government transparency. All information contained 


MN] K N OX herein has been approved for release by the NPS Public Affairs Officer. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 


http://www.nps.edu/library rey Canenies  oe http://www.nps.edu/ 








NPS-62-90-009 


AD~A232 733 


NAVAL POSTGRADUATE SCHOOL 


Monterey , California 


DTIC” 


ELECTE 
MAR 12 1991. 


THESIS 





AN ANALYSIS OF MLAYER: A MULTILAYER 
TROPOSPHERIC PROPAGATION PROGRAM 


by 


Lean-Weng Yeoh 
June 1990 


Thesis Advisor: Hung-Mou Lee 





Approved for public release; distribution is unlimited 
The Chief of Naval Operations (OP-03B) 


Pentagon 
Washington, 


Prepared for: 
DC 20350 








ee re ee 


NAVAL POSTGRADUATE SCIIOOL 
Monterey, California 





Rear Admiral Ralph W. West, JR. Harrison Shull 


Superintendent Provost 


This report was prepared in conjunction with research conducted for the Chief of 
Naval Operations (OP-3B), Pentagon, Washington, DC 20350 and funded by the Naval 


Postgraduate School. 
Reproduction of all cr part of this report is authorized. 


Released bv: 





Dean of Faculty and Graduate Studies 










THIS DOCUMENT IS_ BEST 
QUALITY AVAILABLE. THE COPY 
FURNISHED TO DTIC CONTAINED 
A SIGNIFICANT NUMBER OF 
PAGES WHICH DO NOT 
REPRODUCE  LEGIBLY. 








Unclassified 
secur:ty classification of this paze 

REPORT DOCUMENTATION PAGE 
ta Report Security Classification Unclassified 1b Restrictive Markings 


2a Security Classification Authorny : 3 Distribution Availability of Report 


2b Deciassification Downgrading Schedule Approved for public release; distribution is unlimited. 


4 Performing Organizauon Report Number(s) § Monitoring Organization Report Number(s) 


oa Name of Performing Organization 7a Name of Monitoring Organization 

Naval Postgraduate School (if applicable EC ‘| Office of the Chief of Naval Operations (OP-03B) 

te Address (cim, stare, and ZIP code) 7b Address (cin, state, and ZIP code) 

Monterey. CA 93943-5000 Pentagon. Washington. D.C. 20350 

va Name of Funding Sponsoring Organization | $b Office Symbol 9 Procurement Instrument Identification Number 

Naval Postgraduate School (if apriicadles EC O&MN. Direct Funding 

ke Address (cits. Stale. and ZIP code} 10 Source of Funding Numbers 

Monterey. CA 93943-5000 Work Unit Accession No 
1) Titles inetude security classigeations AN ANALYSIS OF MLAYER : A MULTILAYER TROPOSPHERIC PROPAGATION 
PROGRAM. (Unclassified) 


12 Personal Avthoris) | ean-Weng Yeoh 


12a Type of Report 13b Time Covered 14 Date of Report (year, month, day) 15 Page Count 
Master’s Engineer’s Thesis From To June 199) 158 


In Supplementary Notation The views expressed in this thesis are those of the author and do not reflect the official policy or po- 
sition of the Department of Defense or the U.S. Government. 
17 Cosat: Coaes 18 Subject Terms scontinue on reverse if necessary and identlyy by block number) 

Tropospheric duct propagation, computer program documentation. 





19 Abstract (continue on reverse if necessary Gnd identity by block nusiber) 

MLAYER. a computer program. was developed by the Naval Ocean Systems Center (NOSC) for calculating the signal levels 
of electromagnetic waves propagating in a multilayer tropospheric waveguide environment over seawater. The program is an 
extension of the NWVG which is a trilinear ducting program. Modifications of the XWVG were carried out to handle mul- 
tilaver tropospheric ducts. 

A number of modifications and improvements on the program made over the past several years were not documented. A 
detailed documentation of MLAYER was also not available. The objective of this study is to develop a technical doc- 
umentation for MLAYER using the program as baseline. The study aims to put together the theoretical formulations (spe- 
cific to MLAYER) into a complete self-contained document. This is to facilitate potential users with better appreciation of 
the capabilities, imitations, approximations and assumptions used in the mathematical modelling techniques. As far as pos- 
sible. the same terminologies and functional variables used by Baumgartner (in the XWVG development) and by Pappert (in 
the MLAYER development) are adopted to enable one to relate this document to the program. Step-by-step derivation of 
certain equations was carried out and checked for compatibility with the algorithm in the program. An in-depth scrutiny of 
each program element was also conducted and a description for each is provided. As a result of a detailed analysis of the 
respective algorithm in the program. the documentation for the evaluation of the modal function was eventually prepared. 
Additional materials were gathered from technical reports and papers to supplement the development of this document, 

The MLAYER supporting programs (Microsoft program maintenance utility “makefiles”) were modified to enable the pro- 
gram to run on Microsoft FORTRAN version 5.9. MLAYER was tested and ran successfully on Microsoft FORTRAN 
version 5.0. and C compilers version 5.0. 


20 Distribution Availabihty of Abstract 21 Abstract Security Classification 
WW unclassified unhmned 0] same as report D DTIC users Unclassified 


22a Name of Responsible Individual 22b Telephone cinciude Arca code) ze Ole Symbal 
| Hung-Mou Lee (408) 646-2846 Velh 


DD FORM 1473.34 MAR 33 APR edition may be used until exhausted security classification af this page 


All other editions are obsolete 
Unclassified 












Approved for public release; distribution is unlimited. 


An Analysis of MLAYER: A Multilay er Tropospheric 
Propagation Program. 


. by 
Lean-Weng Yeoh 
Ministry of Defence, Singapore 
B.ENG.(Hons), National University of Singapore, 1983 
MSc(EE), National University of Singapore, 1987 


Submitted in partial fulfillment of the 
requirements for the degrees of 


MASTER OF SCIENCE IN ELECTRICAL ENGINEEP!ING 
and 
ELECTRICAL ENGINEER 


from the 


NAVAL POSTGRADUATE SCHOOL 
June 1990 


Author: a 
Lean-Weng Ycedh 


Approved by: 





Hung-Mou Lee, Thesis Advisor 





Richard W. Adler, Second Reader 





John P. Powers, Chairman, 
artment of Electrical and Computer Enginecring 







Dean of Faculty and Gra@uate Studies 





ABSTRACT 


MLAYER, a computer program, was developed by the Naval Ocean Systems 
| Center (NOSC) for calculating the signal levels of electromagnetic waves propagating in 
. a multilayer tropospheric waveguide environment over seawater. The program is an ex- 
| tension of the XWVG which is a trilinear ducting program. Modifications of the KXWVG 
| were carried out to handle multilayer tropospheric ducts. . 

A number of modifications and improvements on the program made over the past 
several vears were not documented. A detailed documentation of MLAYER was also 
not available. The objective of this study is to develop a technical documentation for 
MLAYER using the program as baseline. The study aims to put together the theoretical 
formulations (specific to MLAYER) into a complete self-contained document. This is 
to facilitate potential users with better appreciation of the capabilities, limitations. ap- 
proximations and assumptions used in the mathematical modelling techniques. As far 
as possible, the same terminologies and functional variables used by Baumgartner (‘n the 
XWVG development) and by Pappert (in the MLAYER development) are adopted to 
enable one to relate this documient to the program. Step-by-step derivation of certain 
equations was carried out and checked for compatibility with the algorithm in the pro- 
gram. An in-depth scrutiny of each program element was also conducted and a de- 

: scription for each is provided. As a result of a detailed analysis of the respective 
, algorithm in the program, the documentation for the evaluation of the modal function 
was eventually prepared. Additional materials were gathered from technical reports and 
papers to supplement the development of this document. 

The MLAYER supporting programs (Microsoft program maintenance utility 
“makefiles”) were modified to enable the program to run on Microsoft FORTRAN ver- 
sion 5.0. MLAYER was tested and ran successfully on Microsoft FORTRAN version 


5.0. and C compilers version 5.0. 
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I. INTRODUCTION 


Propagation of radio waves through the troposphere is strongly affected by the 
complex nature of this propagation medium. The effects on radio wave propagation are 
due to changes in the refractive index with height, absorption by the atmosphere and its 
constituent gases, attenuation and scattering due to climatic conditions. 

Although weather conditions could result in attenuation of radio wave transmission 
through the absorption of energy ov gases, water vapour or scattering by rain, certain 
meteorological conditions can lead to an enhancement of signal strength. This peculiar 
phenomenon known as ducting is due to the inversions of the refractive index profile. 
Propagation in a tropospheric duct can result in an unusual increase in range. 

To predict the anomalous propagation effects, knowledge of the characteristics of 
the refractive index and its effects on propagated fields is required. In the mathematical 
formulation of the model. the concept of modified index of refraction is used to describe 
Wave propagation over a Mat carth. 

In a uniform atmosphere with constant index of refraction, electromagnetic waves 
travel in straight lines over the curved surface of the earth. This situation is unaltered if 
the earth is considered as flat and the rays curved, provided the relative curvature be- 
tween the earth and the rays is preserved. To account for the upward bending of rays 
over the flat earth, the constant index of refraction is transformed to a modified index 
which increases with height and is given by 

mi)= nt , (1) 
where n is the index of refraction, z is the height above the surface of the earth and a is 
the radius of the earth. Equation 1 is still valid in a stratified atmosphere where the index 
of refraction is a function of height. In this case, n is replaced by n(z) in Equation 1. 

The modified index concept allows one to transform a spherically stratified refrac- 
tive structure into a planar laver by preserving the relative curvature between the normal 
of the radio-wavefront and the earth’s surface. Therefore wave propagation over a 
spherical earth can be reduced to that over a flat earth if the index of refraction is re 
placed by the modified index of refraction, For za < 1, the angle of inclination ofa ray 

















at any point in space and the distance (optical length) along the ray path will remain 
invariant in the transformation. 

The flat-carth approximation is adopted in MLAYER [Ref. 1 and 2}. This ap- 
proximation considerably simplifies the mathematical solution of the problem. The ac- 
curacy of the flat-earth approximation was investigated by Pekeris [Ref 3]. It was found 
that for a standard atmosphere, the fractional error in the height gain function is ap- 
proximately 


e= SH, (2) 


Equation 2 shows that the flat-earth approximation will break down at high frequency 
and at great height above the carth’s surface. For example, a 10 GHz signal at a height 
of 2000 m will give a fractional error of approximately 114%. The work of Pckcris also 
indicates that the flat-earth approximation is valid to within two percent for ranges up 
to about one-half of the radius of the earth [Ref 3]. 

As typical valucs of m(z) differ from unity by about 300 parts per million, it is 
sometimes more convenient to introduce the modified refractivity to describe the phe- 
nomena involved in the formation of tropospheric ducts. The modificd refractivity is 
defined as 


M(z) = [m(z) -— 1] x 10° . (3) 


In this program, the modified refractivitv profile of the troposphere is modelled as a se- 
ries of continous piecewise linear segments (see Figure 1). 

The lincar approximation simplifies mathematical manipulations to a large extent. 
In general, any complicated profile can be approximated by several lincar segments and 
the accuracy of the modelling can be improved by increasing the number of linear scg- 
menis and reducing the thickness of each layer. 

If M,(z) is the value of the modified refractivity in the /* laver, then for a linear 
modified refractivity profile, 

M,(z) = M(z) + a (2-2) 3 
: (4) 
2S 2 S ly). 


‘rom Equation 3, the square of the modified index of refraction is 


tw 
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Figure 1. Multiple piecewise linear modified refractivity profile 





me (2) = 1+ 2x 107M (2) — jn , . 


Where 9, is introduced to account for attenuation of the signal due to the absorption by 
atmospheric gases. The derivative of 17? (z) is given by 
fos -6 @M{z) ay 
——{m (z)} = 2x 10 ——— -j— 
dz (ony (2)} dz I Gz (6) 


= a 


Equations 4, 5 and 6 will be used in the mathematical model to represent the character- 
istics of the propagation medium. 

MLAYER also incorporates in the solution of the modal function, the effects of 
surface roughness, atmospheric absorption and the variation of the complex index of 
refraction of seawater as a function of temperature, salinity and frequency. 

The program uses the Shellman-Morfitt complex root-finding routine to locate the 
modes that propagate in the tropospheric ducts. This routine attempts to find all wav- 
eguide modes with attenuation rates below a user-specified value. 

Overflow problems commonly encountered in the evaluation of the modal function 
are avoided by performing numerical computations with extended complex arithmetic. 











Il. MATHEMATICAL FORMULATION OF THE MODEL 
A. HORIZONTALLY POLARIZED WAVE PROPAGATION 


In formulating the problem, the flat-earth approximation with modified index of 
refraction is adopted. The cylindrical coordinate system, (r, ¢, z), is used. The vertical 
height, z, is measured from the surface of the flat-earth. 

A horizontally polarized source may be approximated as a radiating magnetic di- 
pole oriented in the z direction and located at r=Q and z = z,. The magnetic Hertz 


Vector associated with the dipole is given bv 
N(r,z) = M(r,2)2 , (7) 
where F1 (7, z) satisfics ine following equaiion: 

VM (r,2) + k’m*(z) M(r,2z) = —4n {5(x) 6(v) d(z -—27)}, (8) 
where k is the wave number in free space, (.) is the Dirac delta function and z; is the 
location of the source. ‘The right hand side of Equation 8 corresponds to a point source 
of strength —4dz located at r=0 and z=z;. The time dependence e is assumed 


throughout this thesis. 


The electric and magnetic field can be expressed in terms of the Hertz vector as 


—_ 


E = — jou¥ x fi (9) 


and 


—_ 


H=V¥xVxT, | (10) 


where p, is the free space permeability. Solution to Equation 8 can be represented as a 


contour integral: 


:  A(r,z) = F | edott? r., 2), (11) 


“nm 





where H@ (rp) is the zeroth-order Hankel function of the second kind and /{p. z) is the 
height gain function. The contour of integration, C, as shown in Figure 2 can be de- 


formed into C, and C,, also shown in Figure 2. 
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Figure 2. Contour of integration of Equation 11 [After Ref. 1} 





The height gain function in the / tropospheric layer, denoted as f(p, 2), satisfies the 
following differential equation: 


a2 2 | 
\S € i o- = | no z) = —2nd(z-27) . (12) 


Cz” 


In the ground, the height gain function, /;(p, z), is given by the solution of the fol- 


a2 2 
Patiala y=0, (13) 


where #, is the complex index of refraction of the ground and is related to the relative 


lowing equation: 


Ss) 


permitivity, ¢,. and conductivity, o,, of the ground as 


- i ons ee, (14) 
Mags ed eas 2 


The height gain function, f{(p, 2), must represent an outgoing exponentially decaying 


wave as 2 — —o, i.¢., 


f,(p.z) = Ae” , (15) 
where 
ym kala = (16) 
and 
pee. (17) 


For the wave to attenuate into the ground, the branch of the square root in Equation 
16 should be chosen such that 


Infy) < 0. (18) 


If the ground is not dissipative, the branch point and the integration contour satisfting 
Equation 18 is shown in Figure 3. 








Continuity of the tangential components of E and // at each laver boundary re- 
quires that 


fy(0.0) = f,(0,0) , (19) 
E40. = LKG.0} , (20) 
Si (Ps Zin) = Sar (Or 241) 3 P= ltol—-1 (21) 
and 
Eile. tad) = Vin rz) i im dtol-t ey 
Define 
a= [EP bt + ue-2-F] (23) 


and with a change of variable, the solution to Equation 12 can be expressed as a linear 
combination of k,(q¢,) and k,(q,), satisfving 


a2 
E 7 + a| kn(q) = 0 ;n=1,2. (24) 
ce 


This is Stokes’ differential equation. Tor i # | 
Si(p.z) = B,(p) CA; (ek, (gy) + ka (0) (25) 
where k, (g,.) and &, (g) are proportional to the Airy functions 
kg) ~ Ai(—ge"*?) (26) 
and 
ki(q) ~ dAl—4q) . (27) 


Note that the 4, ‘s and B's are different from those in Reference 2 where the height gain 
function is expressed as B,(k, + 4, k;) . Equation 25 corresponds to that implemented in 
the program. 








B — plane 


Figure 3. Branch of Vn? — pf? 


The constants of proportionality of Equations 26 and 27 are of no consequence as 
they allect only Bip) and A(p) of Equation 25 but not f(p, z). The height gain functions, 
J(p,z), used in the computations of the mode sum and path loss are normalized as in 
[:quation 77. The solution of Stokes’ equation is discussed i Chapter IV. 


For the highest layer, i= 1, the height gain function is given by 
Si (0.2) = Brlp)iy (qn) (28) 


where /,(q,) is the modified Ilankel function of order one-third and is proprotional to 


the Airy function as 


hy (gq) ~ Ai(-qe?7) . (29) 





As in Equations 26 and 27, the constant of proportionality of Equation 29 is diso omit- 
ted in the program. 

Lyuauon 25 represents the solution of f(p, 2) for every laver except the one in which 
the source is located. In the laver where the source is located. the height gain function 


is given by 
f(p.2) = Bip) Ly (pe) &y (a + Ala + fp.) , (30) 


Where f(y. @) is the particular solution of the inhomogenous Equation 12. Therefore 


the complete solutions of f(p.2) are given by 
Ap sl = Bip) LA (ky (qd + ke (ad + SO gd) day 5 f= tol I (31) 
and 
fips = Bilpins (qh) + fp. gd bay 3 = 1, (32) 


where 
: l i= NM ai 
CaS | QO GH ° (33) 


It is assumed that the source is located in the 3/* Javer. 
Continuity of flip, 2) and —f(p, 2) at the ground level requires that 
€ 


-° 


le = Bulb Ay (Qi) + Aslai)) + fr Ma) Ory (34) 


and 


a 
al 


. ’ ’ k 3 5 a F 
Jr dy = ByLA A's (gi) + ’2 (4.)) a pea + fo. Mrs - — ¢ 


tae 

, 
Al 

— 


- 


From Equations 34 and 35, 
(fk . 


(a) Part (ng) kL) Beh Bi (Fo ok stawd arhban 
(36) 


= ifn, 914) Os - ae FX. Wrvdowy . 


Simuharly, at laver boundaries, 1 = 1 to [-2, 


10 





Be LA; Aylgeiar) + Antgisent) + SU digs) Oi 
= BC Ay dear ied) + Ads and 
$f dries) Signa 
and 


373 


’ ‘ Z A ina 
BCA ie) +2 Gad] E | i 


cr (5) 7 
+ ae Lf " (p. Tis) On 





> ' ae k 253 
= Bie Cli Adar ign) + ANG ae | ae | Hit} 
es eee , 
eer (P. Giet,i41)] Ota. 
From [Equation 37, 
Ai(Giiay) Bedi + ANG pie VP 
—Kildeay te Bin Aer ~ Aaldigt, 1) Bia 
_ rs) és a ae 
=f (Gist ied) Cet TL UG, te) On 
From Equation 38, 


_\2 k\28 
Aq: wd ee) 4; BA; + Ry q;, uo #) * a, B, 


{ 





4 k 2/3 
— K's (digs, al ) 4) Bigg ligy 


O44 


k 2/3 
— kyxgia), ma ( ) O41 Bin 


Hitt 
a {5) 5 
= ae PP. digs, 41) 5 Oia 


C ‘ 
a ps Uf (p, qj, it1)] On. 


(39) 


(40) 


Rig Lap Aylgien) + ld I tL OU a D8 aay 


(41) 
= Bh land + fOr, qineiny 
and 
. 7 7 k 23 
Brey CAA (qian) + Kd DI as a es 
fa 2 
ae C6. 91-1.) 10 7-1,49 (42) 
' k t C S ee 
= Bh staun( x ) °ap + = dak CO) 
From Equation 41, 
Aid) Bey lye + AQ) Brey — hl Br 
{s) : (5) ¢ (43) 
=f (peg ders mL (Pe dar DOr yay. 
From Equation 42, 
<5 Kk \23 3 kee Noug ; 
K'i(dqypy cae Oy) Byy Ayey + hal qeey yp) ay) J OMlat Byy 
; kK \2/3 
— Wx) a) B, (44) 


= = Le tes 9119144 ~ = tp, GH-ypeay 
In the above equations, q., is defined as 
Uae = (z= 2) . (45) 
Equations 36, 39, 40, 43 and 44 can be written in matrix notation as: 


Ag =4, (46) 


~ PS 





where 


As, As, Az; 0 
A;, As; Ay, O 
A, Ass Ags Ass 0) 








OQ Aspages Ayesgra Aseae as Ayesci2 0 
OQ Agsars Asara Ansara Arsana 
0 Aspsa-3 Azer: Areas 
v Arraga3s Axeisrez Azra 
(47) 








Bats 
B, 
fos ; (48) 
Bey Alyy 
By 
B, 
and 
Wy 
Wy 
y=: : (49) 
Wy 
The elements of the matrix, A, are 
kK \28 4,; i - 
Ay = ee OA Gy) ~TA Ia) (0) 
be (NI ie noes . 
A\2 = Oy ak 91.1) —Jrk2(93,1) ’ (31) 
A>) = ki(q)2) » (52) 
Ay2 = A2(92) 5 (53) 
A33 = — ky{q2 9) 5 (54) 
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Ars = —Arlg2) . (35) 


- 


. ko \23 5 
Ay, = & wna 4 ) Hy (56) 
Ay = Mana = ya ; (57) 
Ay; = - Klaza( ype Oy (58) 
Aaa = kaa $Y? 2s (59) 
Aspazies = ANG) + (60) 
Ayagrn = A(Qei) 5 (61) 
Aypaze1 = — Ald) 5 (62) 

, k 2/3 
Ayi-121-3 = k (ad) PEeh (02) 


' kK \2;3 
Ax-121-2 =k 91.9( GES ) ae 








, kK \23 F 
Aycan OA Aaai( 45 ys (65) 


The elements of the column vector, w , are 


~ ls) , a 5) i 
y= pel. AM a ~ FPA ist + (66) 
ea fp, 922)02\1 — fp, Nv (67) 
C ) : @ seh . dpe 
ay 1, 4220182 - ZU og aiia « (68) 
Wi = fp. q1.)8 114 = {Cp T1817 (69) 
and 

C : 7 c ‘) 3 - 

v= ao le Hi ~ 3 “p, Fr DIO nay (70) 


To determine the modes that propagate in the tropospheric ducts, the roots of 


|All =o (71) 


have to be located. This is discussed in Chapter III. The symbol. || A ||, denotes the de- 
terminant of A. 


From Equation 11, the solution of Equation 8 in the / laver is 
| ry QO) 7 Z > 
Mrz) = +] pdpHy” (rp) filr,z) (72) 
ad 
From the theory of residues, Equation 72 can be expressed as 


16 








4 


BE RO gee. oe EY Supe aie Oinrc der <8 
ZPitm) +5 eet (rp) fi (p. 2) 


> 
< 


Nurwc) = - 





+ +| papll, @) (rp) fip.z) » 
¢ 


- . 
5 
2 


where 6, (m) are the residues of p//, (rp) f (p, z) evaluated at the poles of f(a, z) in the 
fourth quadrant. The contour C, is the quarter-circle of infinite radius and C, is the 
contour enclosing the branch cut of the radical of Equation 16. Since /73' (rp) decays 


exponentially in the lower half complex p-plane, 


| pdplly i (AF (p.z) = 0. (74) 


Cy 


For small height. z. the integral over contour C, represents the contribution by the sur- 
face wave to the total field. Bevond the horizon, this field is small compared to the total 


field and can be neglected. Therefore Equation 73 becomes 


Mu.z) = - nj) by Qn). (75) 
Mn 
Where , (m) is given by 
bm) = Hy Pm) Bm (2) 8m Er) (76) 


and g,, (z) is the normalized height gain function given by [Ref. 2 ] 


Jnl q=)] 


8n(2) = ’ (77) 





— Sm CHO fds” 
ads ee 
I Sin [9( )] “Pm (J dp io Pm 


where g(z) = q(z) for z, < z < z,.;. 
According to the program source code, the integral, // = fe *(q(z)]dz, is evalu- 
ated by considering that f.[q(z)] is a solution of Stokes’s equation, 
hig 


=a din = O . (78) 
dy 











This leads to 


_d_ r 2 dh, _ 2 ie dis. hn a hin 
dy: +4 dg | =fin + Ly Sin dy Te dy : dy" 
din [ttn (79) 
= fn” +2 | OG 
q dq 


y 


=Sm> 
Note that Equations 78 and 79 imply that the derivative of f, is continous at z;. 
is. f(g) does not include the contribution from the last term of Equations 31 and 32 


The integral of f, *[g(z)] in the @ Javer is given by 


That 


Ppa ) ] Q-) 7. 
ae Cy(< J d= q'; hy; (q;) ag; 
=; t"¢@., 


} (SO) 
df, (qi) |2 | Gant 
= yy i (;) +| oe | a ’ 
where 
tq: 
20% (81) 


From Equation (S80), it cun be seen that the contribution of A/, to // atz = z, (<1) 


ss cle) dfn (gd 2] 
Al, = D gilin (gi + | a a 


l 


is given by 
(82 


| 
pms aaa dq; CS nt + ie gigt CSmt qi- Di 


df. 
By making use of the fact that ae is continous at the laver boundaries, 


ns hn (Ga) = ae Lon (4) 





which imphes 


4 
tei ae Tye Ldn | G42 = 12g, Ln G9] ($4) 
and 
f GL Ge es 2 
‘oes rere Les , = i § 
t 7 ~ Em | a we = ( Gin! } | dq; Ln: gM) (85) 


- Substituting Equation $5 into Equation §2. A/, becomes 





i A ($6) 
| ty") | d 2 
+| --—- = Aga 
q'; ia alae ing } 
For1=1}, 
Al, = =e aad n=] Ua tol : (87) 
GP Gag at 
The integral, //, is 
t 
i= > Ay, (88) 
er 
The second term of the denominator of Equation 77 can be written as 
jee CL foe ake fo TOT fds dy 
“Via i{ dp Vp. — 2 pix ak dg, dp yee 
fn THO a k \23 2 Pm 
ais tg ¢ (+) t (89) 
ol dq, " i 


_ fn La] ( ke Ne dl 
ae kK? a dq, 


For cases of interest.|p,7{ > 1. The Hankel function in Equation 76 can be replaced 


by its asymptotic approximation: 
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2 2 hal : 7 ; 
Ho (pr) ~ Sao el Sy ae (90) 
‘stn | 


Therefore Equation 75 becomes 


2 2 : ae 
Nr. z)=— yy (a> Ni exp ( Thm a iar ) m=) Em (27) 


a eT SNR ON SB ee ee dig 
re Tr e gk &m (z) &m(27) e (Y1) 
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To allow for spherical spreading. the r-!? term in Equation 91 can be replaced by 


. r ' : ‘ ; 
Ca sin (7 )J??.a being the earth’s radius. Equation 91 becomes 


pay 1/2 i -fass + =132 Hpi pt “4 

Mrz) = Edie epee, e0 Pm’ &m (2) &m (zy)e . (92) 
asin(—) Af 

d m 

In the cvlindrical coordinate system. (r. @. z ), the transverse clectric field component, 

E.. is given by 


a 


be = — wy a ; (93) 


Differentiating Equation 92 with respect to r gives 


fale Be cos (77) or ames cree 
Cr ad 2a vsin'( + ) asin(—) m 
(94) 
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Neglecting higher order terms of = . Equation 94 becomes 
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The free-space solution to Equation § is given by 











jr 743 —jer 
pI ee gu 
ps === = : ,r>e:s. (96) 
\ reer 
Thus, 
€ IT S ik ~j); —ikr 
— = ie (9°) 
cr 4 = 


ae : 7 ae : : Seo 
and similarly. the higher prder term of — in Equation 97 is neglected, giving 


cn 1s Jk 


cr a E 


ee (98) 


Normalizing the field to the free space broadside field, the mode sum becomes 


dey? Vee x ee oor 
es = a ' ~ - a Ibs 
a ai Pm Em (=) &m ( 7) oF : (99) 
“OFS ka sin{ =) Mm 


The coherent mode sum in decibels is given by 





ECMS = 10 log 


(100) 
f 1/2 TP mi? ‘ 
= 10 logy; ———— 7pm 8m (=) &m (=z) ’ 
l k°a ae J 


The power sum or incoherent mode sum (EIMS) signal can also be calculated if the 


phase of each term is ignored: 


=p pt | 2) 


pe f Qrr? Ss 1)2 a 
ELMS = 10 106 Min Cele) Sakae ae (101) 
eae se j 







The path lossses are 


Coherent path loss = 32.454 20 logr+ 20 log f-ECMS (102) 








I 


32.48 420 logr +20 logf—EIMS . (103) 





Incoherent path loss 






Where ris the range in kilometers and fis the frequency in MIIz. 






B. VERTICALLY POLARIZED WAVE PROPAGATION 


A vertically polarized wave mav be treated as due to a radiating electric dipole on- 






ented in the z direction and located @ r=OQ ands = c.. 


The electric field and magnetic field can he approximated by [Ref 1] 






H = jtgia¥ x il (104) 





where [7 is the electric [lert7 vector. 


Calculation of the electric [lertz vector is identical to that of the magnet Hertz 





vector except for the following modifications: 






1. The height gain function satisfies the following diflerential equation: 







2. Continuity condition at i= 1 becomes 






nehlp, 0) = NOW (p, 0; . (107) 


This reads to 








134, ; nO) , 
\P aA Ndi mie a Aa) (108) 





and 


' my iO) ~ (p) ~ c (p) z 
le > ie — ops ali — IT enya (110) 


Hy 


The transverse magnet field component, //. , is given by 


a 


IT, = = JOE call . (1 ] 1) 


oF 


Prom Fquations 94 to 9S, it can be seen that 








Il, Es 
fps . Lets 
(112) 
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Pad Ve iz spat 
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It follows from Equation 112 that ECMS, EIMS. coherent path loss and incoherent path 


loss are given by Equations 100, 101, 102 and 103 respectively. 


C. DETERMINATION OF THE HEIGHT GAIN FUNCTION 


The height gain function, /(. 2), 1s given by 
inc) = BCAA Mg) + hotg)l oo: i= ttol—1 (113) 
or 
ii(p.2) = Bh(q) viel . (114) 


The coefficients, -f, and 2, are functions of p as in Equation 25. 

Determination of the values of 4, and B, is dependent on the evaluation of the in- 
tegral (or summation) of Equation $8. It was found that if the attenuation rate of a 
particular mode is small (< 0.1 dB km), the integration Equation $8 is carried out from 
the top laver to the first laver with ihe height gain function at the «*p laver set to a 


certain finite value. This will prevent the integral from blowing up. 





On the other hand, if the rate of attenuation is large (> 0.1 dB km), the integration 


is performed from the first laver to the top layer with the height gain function at the first 


laver set to a certuin finite value. This will ensure that the value of integral will not be- 


come too small. 


Therefore, for attenuation rates greater than 0.1 dB km, 4, ‘s and B, ’s are deter- 


mined successively starting from 4, and B, up to 4, and B,. Ati = [andz = 0, 


LS. ee Q) = Si(p. r=) 


and 
LAP. Meco = Gr CA. ae 
Thus, 
Ay = By CAA a) + Ala) 
and 
; a pa, of A \2;3 
Fide = By TAA ga) + 4 9.0 # ) oy 


Dividing Equation 118 by Equation 117 gives 


kK \273 4, ; 

( oy )ak AM) — FAG 1) 
A4,=eOororoorr 
ET Ke NDR a 
JAM) > eye ok (41,1) 


Q'y k'ldy yy) — Frag) 
WAC) — TAG) 


where 


» _ 4h = { -& \23, 
ee dz a 4 


In the program, /((p, 2,) is set to 1. This gives 


PCRS IET ns eee 
ACA (Gy 1) + Anta 1) 


(115) 


(116) 


(117) 


(11S) 


(119) 


(120) 


For 1=2 to I-l. f, ’s and B ‘’s are determined by matching conditions at everv laver 


boundary. At laver boundary z = 2, 
fAp. 2) = flp. =) 


and 
= lee = Sr Late. -i 


- From Equation 122, 
BCA Agi) + An(gis)] 
= Bola Agi) + ANG 5 1 
Thus, 


Bi CA As diy )+ klg_y 3) 
AK (Gi: ) + ANG: ) 


B, = 
The coefficient, .f,. is determined from Equation 123 as follows: 
BLA KR ga) +A AG 07 
= Bolte Alger + haere die. 


Substituting Equation 125 into Equation 126 gives 


Ary Aylgiy + ANG) ‘5 
ae [AA (qis) + A olga id’: 


Pein A Gera) +A aged at 


Thus, 
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ny; 
de Ee 
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(122) 


(124) 


(126) 


(127) 


(128) 


where 


nA, = UK Aga ding AG) + AQ: yj J 
(129) 
— Oy Age elie A Giera) + alga) 
and 
d4;, = iy Ayia Ain A'vqira) + K'(Gi-4 3 )J 
. (130) 
—P KG AR Ay(Ginn y+ AG) 
Vor i = | and at the laver boundary. z = z,. 
Si (pet) = fray (21) - (131) 


This gives 


Rive Bri Celney Avie) + Ad] (132 
oo Aldyy) = 


Dor attenuation rates less than or equal to 0.1 dB km, 1, ‘$s and B, ’s are determined 
successively starting from A,and B, down to 4, and B,. 


The height gain function at the top laver is given by 
Sylp.z) = Asta). (133) 


Where the coefficient. B,. is set to one. Matching the boundary condition at z =z, gives 


Boe hylqyy) (134) 
Arey Ayldien yy) + Adlgeas) 
and 
AA; af 
Aj = 7 Ti,’ (135) 
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where 


Ap = YW ap) Aq) 
| (136) 
— fy Milan Kg 1) 
and 
AA = Oey ly R911) 
(137) 


= hai Aa 1) 


Similarly, the coefficients. 4, and B, are obtained by matching boundary conditions at 


z=z,fori= I-2 to 1. 





Thus, 
A; k (4; i ) + Kilq; ; ) 
Fiiiig +1 Tied aaa: +1 itd (138) 
AA (dicey) + AAG 4 ) 
and 
pe 139 
a te te) 
where 
WAG = digs Ao(Giigy igs A Gee ign ) + A294 ar J 
(140) 
— PR dig Clin Avert) + Ae 
and 
AA, = QTR Gi iy Aig Ai CGig egy ) + AaGig iia] 
(141) 


— Figs Apia, Aig A's (Gist ie) t+ Aa(Qigs ver JJ 


D. TEST FOR EVANESCENCE 
The functions, 4,(g) and 4g). in Equation 26 and 27 can also be expressed in terms 


of the modified Hankel functions of order one third as follows (sce Appendix A): 





, : , 4 iy 
= { srry ane ($9") (142) 
= hy (q) 
and 
Ax(g) = 2(12)1" e/7° 4( —9) 
(143) 
=hy(q) — eng) , 
where 
4) : t . + 
hig) = Cea ae ($0°") (144) 
and 


145) 
— { 2 32\123 Gy fe 22.318 ae 
= ( 3 q ) My 15 ( 3 q ) . 


Tor large | q|. A,(qg) and &,(q) can be approximated by the asymptotic expansion of the 
Hankel function (Ref. 4 and Appendix A]. That is, 


ts 5.5 as eo 
kg) ~ (a aH 8 exp (Sg — 12 )| 


2r dz 
at arg(q) < aa 





(146) 


and 


16 Wg (2 §, 
kylg) ~ (12) a? g exp] - ($a 2 7) 


ty 


(147) 
O < arg(q) < dn. 


From Equations 146 and 147, it can be deduced that for large | g| and in the region 


ty 
too) 


I~ d- . 
— < arely) < —, (148) 
Pea 


A(qg) will be exponentially large and &,(qg) exponentially small. When this occurs, 


A(qg) and A.(g) become linearly dependent numerically, ie.. from Equation 143, 
hig) = e*8halqg) . (149) 


This implies that A, (g) and A;(g) are also linearly dependent numerically. Under this 
condition, the matrix of the modal function will become singular (in a numerical sense) 
for all values of g. Pappert and Goodhart defined the fields as evanescent when Equation 
149 holds [Ref 5 and 6]. 

When evanescence occurs, the height gain function is represented by an exponen- 


tually decaving field, Le. 
Silp.=) = Byky (qi) . (150) 


It is necessary to incorporate test conditions in the program to check if the rf field 
has become evanescent. This is achieved by tracking the real part of the exponent of 
k,(q) in Equation 146, 

Since evanescence may onlv occur in the region specified by Equation 14%, it can 


be deduced by simple transformation of Equation 148 that 





: 2. 213 : 5 
“—— < arg Id! < aA (141) 


This shows that one needs only to test for evanscence when the real part of the exponent 


of k, (g) is positive. That is, 


Rel 23 q"?| SO 4 (152) 


As the coefficients of the height gain function are dependent on how the integral of 
Equation 88 is evaluated, it follows that different test conditions exist for upward and 
downward integration. For upward integration. conditions for non-evanscence in the i” 


laver are 
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2 Raa, 42 
Re | $7 (ica) 2] = Ref Sst"? | <7 
= 2. 
for Ref Sil aa | a (153) 
ae 372 
and Ref z/ gu? | > 0 
or 
2 32 2 3/2 
Ref Si ised)” ] Sa for Ref Sita? | < 0 
(154) 
; 
and Ref Silas) | >o0 
or 
dq; 
Re > 0 (153) 
dz 
or 
-i41 < 45 on . (156) 


In principle. the right side of Equation 153 should be as large as possible. However, a 
large quantity will lead to severe cancellation error when extracting 4, (g) (which is a 
small number) from a linear combination of two large numbers, f, (qg) and A; (gq) . 

Equation IS4 is obtained from Equation 153 by setting Ke[ <j (q_.)87] to zero When 
its value is negative. This is true because when the real part of the exponent of A,(q) is 
negative, the condition of evanscence will not occur. 


Ps cae ; ad adn. 
With regards the condition in Equation 155, it can be shown that when —— is small 


(which is usually the case) and Equation 155 is true. then =< 
Re(a > 0. (157) 
This implies 
me? (Za) > (a) (158) 


which leads to 











Reldieg] > Relg (159) 


or 
: ae 3/2 re 6322 . 
Re wane < Re AC ae (160) 
or 
2 Ce ee: 
Ref Zi gua” = Re| Sita.r"* | a8 is (161) 


By comparing Equation 161 with Equation 153, it is clear that Equation 161 is a condi- 
uon for non-evanescence. 

Finally, the condition in Equation 156 was based on the fact that at lower altitude, 
the magnitude of q is small enough such that /,(g) and /, (g) are well behaved in the re- 
gion stated in Equation I4s. 

By the same argument. conditions for non-evanescence in the /* laver for downward 


integration are as follows: 








. Be 2 rie 3/2 
Re a (917-4) ve = Re Tad (g:_) 3) ie < 7 
3 2 2/2 
for Re | $0 (Geel | = 0 {162) 


5 
and Ref Fi tain, ta > 0 


or 
ey >) 4 
3 ; (163) 
2s 3/2 
and Re [WG-11 pe ee 
or 
Ay 
. Re| —= Se | < 0 es) 








or 
z= < dim. (165) 


E. EVALUATION OF THE MODAL FUNCTION 
When evaluating the relative ficld strength of the electromagnetic waves using 


Equations 100 and 101, the eigenvalues, p = p, , for which 


|All = 0 (166) 


are required. The determinant, || A ||, is known as the modal function. Each eigenvalue 
represents a different electromagnetic mode that propagates in the multilaver tropos- 
pheric waveguide environment. 

The roots of Equation 166 are found by the Shellman-Morfitt complex root-finding 
routine to be discussed in Chapter II]. Instead of searching for solutions of Equation 
166 in the p-space, it is more convenient to conduct the search in the g,,-plane. The 


variable, p, is related to g,, by 


k 2; 2 Pn \2 5 
1, = (a) | my - (4) : (167) 


In searching for the roots of the modal function, the determinant must be evaluated 
many times. Therefore an efficient method for accomplishing this is required. In partic- 
ular, the Laplace’s expansion is used to obtain the determinant. The Laplace’s expansion 
is based on the following theorem [Ref. 7]. 

Theorem: Laplace’s Theorem states that if one selects anv r rows of a de- 
terminant, |.{|. forms all possible r-rowed minors from these r rows, multi- 
plies each of these minors by its co-factor and then adds the results. one 
obtains | 4 |. 

The method also takes advantage of the presence of zero elements in the determi- 
nant to reduce numerical computations. 

The understanding of this technique can be best achieved with an example. Con- 


sider a trilinear refractivity profile with the following modal function: 


es) 
ro 
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[PA ly = 751 7 ae 3 : (168) 
Q 0 300 gy Gs 
0 0 Qs3 0 sy ss 


Let the last two rows be the r rows from which minors are to be formed. Then applving 


Laplace's Theorem, the modal function is given by 


4, 42 @ 


AY; = [ea ee 33 
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ayy ay 0 
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ee] fa 2 23 
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Ifa four-laver piecewise linear refractivity profile is used. the modal function will be 


a4 @ OO 6 0 oo of 
of 2 az a3 G4 (Q) Q) Q 
Cay a> 433 44 Q Q 0 
Al, =] 0 0 43 ay4 ays U4y 0 . (170) 
0 0 33 a3 ass Usg Q 
0 Q Q) 0) a5 Qov5 a> 
Q) 0 0 0 Ors. ~ sey. <4ae 


Notice that Equation 170 is built from Equation 168 by adding eight non-zero entries. 
NamMelv. Ay. Gee Chee Gege Ayre Cees G, and a--. Therefore, the results from Equation 168 


can be used to compute Equation 170, ie., 











a % 0 O DO 
>; G5. Oy a 
Gg 57 Oo ee ee | y es tae = 
fA, = WAN; Gout = 3, 32 433° Bd Q) ie ae: (171) 
0 Q 33° 44 Aye 


In Equation 171, the determinant of the 5 x 5 matrix is also evaluated using the 
Laplace’s expansion. 

Without loss of gereralitv, the evaluation of modal functions of multilayer piece- 
wise linear refractivity profile is performed bv systematically increasing the number of 
rows and columns of A and applying the Laplacc’s expansion. 

The derivative of the modal function. || A |], with respect to g,, is also required by the 
Shellman-Morfitt root finding routine. This is accomplished by expanding || A | using the 


Laplace’s expansion and performing differentiation using the chain-rule. 











F. EFFECTS OF SURFACE ROUGHNESS 


The materials in this section are adapted from Chapter HHT of NWVG [Ref. I]. Tuas 


included in this document for completeness. 
Surface roughness could have significant effect on the signal levels at large ranges, é 
especially when the frequencies of the rf waves are above several gigahertz. 
The effects of surface roughness is modelled by Baumgartner in NWVG [Ref. 1] as 
a lossy infinitesimally thin bottom laver of constant modified refractivity shown in Fig- 
ure 4. 


In the bottom laver.0 < z < z,, the height gain function, f (p, z), can be written 


as 
fole.z) = Ag fel + Re (172) 
where 
L= ky n(O)— fr, (173) 
: pr s 
De ame (174) 
ko 


and R is the plane wave reflection coefficient. 
By matching boundary conditions at cach laver boundary as was done in Section 


A, it can be seen that all previous results still hold if the following substitution 1s made: 


1—R 
1+R 


— Yo= uu 


sy 
~ 
i 


In general. the bump heights of the rough surface are assumed to be Gaussian dis- 
tributed. In this case, the reflection coefficient, R, can be represented as the product of 
the smooth surface plane wave reflection coefficient and a surface roughness factor [Ref. 


8]. For horizontal polarization, the reflection coefficient, R, is given by 


4 x 2'3 
Ry ex =24°5°( 4b) ; Re > 0 
R= i. u P| As ri M1 (9),1) ; (176) 


where 6 is the root mean square surface bump heights and Ry is the Tresnel reflection 


coeflicient for a horizontally polarized wave given by 
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Figure 4. 





Effects of surface roughness modelled as layer 0 [From Ref. 1] 
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For vertical polarization, R is given by 
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a ea 23 5 
' ( Reeve) -24°5°( 3 } te > Re(qy 3) = 
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(17S) 
l Ry 2 Relgy yl) sO 
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where &..is the Fresnel reflection coeflicient for a vertically polarized wave given by 


ne on (O) = py? = ne (0) [1 — rel * _ 
Ky = Fa a A a eee i (179) 
LN ae ae ai OL a fe 


For horizontal polarization. ;, can be expressed as follows {Ref. 1]: 


alg, ,) 
; b (dy) 





sn = (180) 


where 
ald) = Ga a aa tanh (+) + k + (=. qr i . (181) 
b(q)4) = 14 ae (gay? tanh (S| T4 Cou al . (82) 
T=.» 1” (0) (183) 
and 
b= 275 (BY) ae (184) 


For small values of | q,, |, @(9,,) and 6(q,,) can be expanded in series as 
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The derivative of ;., with respect to g,, can be obtained bv dit’erentiating Equation 


ISO. Le. 
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Similarly, a and a are obtained from Equation I8] and 182 respectively: 
qi ain 
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For small values of |g, |.. 
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The corresponding results for vertical polarization are: 
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For small values of |g.) |. a(g,,) and D(g,1) becomes 


39 

















~ op AES V2) p252¢ ZY 273 a 
ag, 4) = Se) qh) or fs a 6 Geert ok Cee) q) 
Jee GAO AE Oe, 
+ 15 k () ( k ) Gy ne eeee i 
nO) 973 1 1 4.2/3 
Se tap hye a, (196 
! 1 3/3 2 ] 1 \2 3 
Sea ge qi a (FV a 
Ss 8/3 7 M1 10/3 5 l 
= ) + : lees 
ps Pee me ee S 
and 
m(O) 4, y y 
b(q, 1) = | —_— : peer ae i ‘. 7) arate + 5 ai 
6 ga o 
~| yas SS 6" | (Sh) 3 Yh) 
3 ‘ k 
eee eg eye (197) 
6T oT ee se = 
2 10 -10 1 HF é oy 3 4 
emcees Lo =a 
+[ 6 Page eg |pta 


The derivative of a(q,,) and bag) are Obtained from Equation 196 and 197 respec- 
uvely: 
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and 
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For small values of | g,,/. == and = becomes 
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G. COMPLEX INDEX OF REFRACTION OF SCA WATER 

The materials in this section are adapted from Chapter IV of XWVG [Ref. 1]. 

The complex index of refraction is one of the physical parameters required for 
evaluating the modal function and height gain functions. The complex index of refrac- 
tion of sea Water, n,. is a function of frequency, temperature and salinity. It is related to 
the eflective relative dielectric function, e,,,. and the effective conductivity, 


a... through 





= et x”)? 
Ng -- Cot J Cy : (202) 


where €, 1s the permittivity of free space and w is the frequency at which #, is evaluated. 
The effective relative dielectric function and the effective conductivity are related to the 


complex dielectric function, e, and the ionic conductivity, o. through 


e= eg — je", (203) 
fo = € (204) 

and 
Gey = T+ Wy ae (205) 











In Equation 205, @ is assumed to be real. This is a good approximation for frequencies 
much less than 10? Gllz [Ref 1]. 

The theory of an ideal polar dielectric in an alternating electromagnetic ficld was 
first formulated by Debve [Ref. 9]. In the theory, the complex dielectric function is ex- 


pressed as a function of a single relaxation time, t, as follows: 


e=e—je’ , (206) 
Ce ek 
es — ‘ (207) 
L+ oot 


and 


(6, — € det 
f= OO (208) 


l+oor 
where ¢,, is the part of the dielectric function due to the atomic and dielectric polariza- 
tion, €, is the static dielectric function and w is the frequency of the electromagnetic 
wave. Equations 206 to 208 represent the fall in the value of the dielectric function from 
e.to é,. This fall is accompanied by a single broad absorption band in the vicinity of 


the characteristic wavelength, 4,. given by 
; | 
fc eS (209) 


The characterisation of the variation with frequency of the complex dielectric function 
of water in terms of a single relaxation time is valid for frequencies up to 300 GHz 
{Ref. 1]. For higher frequencies. more than one relaxation ume may be required. 

In the model of Klein and Swift [Ref 10}. c, is assumed to be 4.9 for both pure 
water and seawater. With this assumption, the static dielectric function for pure water 


is given by 
tsp = 88.045 — O.4IA7 1, — 6.295 x 107 2,7 + 1.078 x 107 1? (210) 


and the static dielectric function of seawater is given by 


ess = ts, a(S, iE) ’ (211) 
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where 


esq = 87.134 — 1.949 x 1971 7, = 1.276 x 107 29E x oF 12) 
and 


a(S. 1.) = 1.0004 1.613 x 107° Sr, — 3.656 x 107 § 


+3, 107? S$? 432% 10° 


In the above equations, 7, is the temperature in degrees Celcius and S$ is the salinity in 


grams of salt per kilogram of seawater. The relaxation time of seawater is given by 


RST): as (214) 


where t, is the relaxation of pure water given by 


t, = 1.768 x 107! — 6.086 x 107 1 L104 x 107 7? = SD x 107 a? (28) 
and 
A(S, 1.) = 1,000 + 2.282 x 107° S 1, — 7.638 x 107 S Oe 
a 216) 
— 7.760 x 107° S74 1.108 x 108 S* 
The ionic conductivity of seawater in siemens per meter is given by 
o(t.. S) = o(25. S) exp(— Ag) . (217) 
where 
A= 2-1 , (218) 
@ = 2.033 x 107? + 1.266 x 107° A + 2.464 x 107° A? 
. (219) 
— S[1.849 x 107° — 2.551 x 107° A + 2.551 x 107° A] 
and 
o (28, S) = S$ [0.182521 — 1.46192 x 107° § + 2.09324 
(220) 


x 107° §? = 1.28205 x 107° S$? ] 
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H. ATMOSPHERIC ABSORPTION 

The materials in this section are adapted from References 11, 12 and 13. 

The gas molecules in the air absorb energy from electromagnetic waves. The ab- 
sorption is due to the quantum mechanical resonances of the gas molecules. These res- 
onances at specific frequencies then broaden into bands as a result of collisons of the 
gas particles. In the lower atmosphere, the density of gas particles is higher. This means 
that collisons are more frequent and the absorption bands are broader. 

In the troposphere (height < 10 km), only the lines bands of oxvgen and water 
vapour are important. The principal resonance of oxvgen is at 60 GHz and that of water 
vapour is at 22.2 Gllz. The extent of absorption by oxvgen molecules is dependent on 
the pressure. p. and temperature. t. In the case of water vapour. absorption 1s also de- 
pendent on the water Vapour density, ». For specific values of p. t, and p, the absorption 
coeflicient. g., .1s computed as [Ref 11] 


an = WIS20{N, —— dBikm_, i221) 


where fis the frequency and NV, is the absorption spectrum. The absorption spectrum, 
NV. . is related to the line spectra of absorption, SI°. the continuum spectra. .V,. and the 
liquid water extinction, .V, . by the following expressions [Ref 11): 


Vyo= SF + NM, + pp. (222) 
The summation of St 1s done over 44 spectral lines of oxygen and 29 spectral lines of 
Water Vapour, [Ref 11]. ie. 
da 


SISF = > a) (i) pO? exp [a () (1-4) F 
i 


i=] 
+ Db (ice? exp[hy)(1- OF , 


Where a,(/), a, (/). 6,(é) and 6, (/) are spectroscopic coefficients derived from experiments 
and spectroscopic parameter compilations. Values of a,(7). a, (4). 64) and L(A are listed 
in Table J. 











The drv air pressure. p. at a height of z above the earth's surface is given by 
[Ref. 12] 


p= 101.325 1 re a kPa. (224) 


The relative inverse temperature, @, is related to the temperature, ¢,, in °K by 


se 
= ao: (225) 
we 


The water vapour partial pressure, e,. in kPa is related to the saturation water vapour 


pressure, e,, and relative humidity, rh. by the following equations: 
=) 
ey = e,{rh)x lo™ , (226) 


mye . 
€, = 0.6105 exp [252 ( ie) — 5.31 log ( 315 )| kPa (227) 


and 
i= t — 273.15°C . (228) 


In Equation 223, F is the shape factor of the broadened resonance lines given by 
[Ref. 11] 


—— (5 Le iia + gg Bone | (229) 
ie OA ona 2 i oe (vot Ni +e 
for oxygen line I to line 38 and 
f* g g 
fe (=;)| a= ‘ — | (230) 
ee aaa aie 4 iigcky) ee 


for oxygen line 39 to 44 and the water vapour lines. 
In Equation 229 and 230, v, is the molecular line centre frequency, g 1s the width 


and ¢ is the overlap interference. 











For oxygen in air. 


g=al(ptl3ee” , (231) 
€ = a, pe" (232) 

and for water vapour in air, 
g = b,(48e + p) Hee (233) 





The line data a,. a, and 4, are listed in Table 1. 


The continuum spectrum, .V., is given by 


Vo = 19 fep® x 107° + 6.28 ppd? x 107 [i + 2.1pe'* x | ppl 234) 


™ 
Pay 
+ 
Le) 
=r 


and 
gy = 0.012 s 238 
ayo . 2p + 1.3 e) 0 A (235) 
The liquid water extinction. .V, . is evaluated from the following expressions: 


gy 
Nee es. priiony. Sans (236) 
(eo + 2)° +(e") 


and 


Nye = 0.55 wf! o* ppm for f > 300GHz , (237) 


Where w is the liquid water concentration in grams per cubic meter. The real and im- 
aginary part of the dielectric function of water, e’ and e”’, are given by Equations 207 and 
208 respectively. It is assumed that ¢, and given in Equations 207 and 208 are given 
by Equations 210 and 215 respectively. 

The absorptivity of tropospheric gases, 7, is related to the atmospheric absorption, 
g.5, in dB km, by the following expression [Ref. 13]: 


ered 
Pers 


Sa = 10° logge . (238) 








Hence, 


pe re. (239) 











Table 1. 


SPECTROSCOPIC COEFTICIENTS FOR OXYGEN (a,— a) AND 
WATER (6, — b,) [FROM REF. 1] 
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Center 
. frequency, vy Strength, a, Temperature Width, a, Interference, a, Temperature Quantum number 

(Gliz) (kHlz/kPa) exponent, a, (Gllz/kPa) (1/kPa) exponent, a, identification, K * 

0" 0.307 E-3 (ppm/kfa) 0.12 (0.056) E-I I* to 37” 
50.4738] 0.94 E-6 0.969 El 0.86 E-2 0.520 E-2 0.179 E+t 37° 
50.98742 0.244 E-S5 0.869 E+I 0.87 E-2 0.550 E-2 0.169 E+l 35° 
51.50311 0.604 E-5 0.774 E+1 0.89 E-2 0.560 E-2 0.177 E41 33° 
$2.02124 0.141 E-4 0.684 E+] 0.92 £-2 0.550 E-2 0.18) E+i 3 
$2.54227 0.308 E-4 0.600 E+I 0.94 E-2 0.569 E-2 0.179 E41 29° 
53.06683 0.637 E-4 0.522 E+% 0.97 E-2 0.528 E-2 0.189 E+] 27° 
53.59570 0.124 E-3 0.448 E+1 0.100 E-I 0.544 E-2 0.183 E+t 257 
$4.12997 0.2265 E-3 0.38! E+ 0.102 E-1 0.480 E-2 0.199 E+] 23° 
$4.67115 0.3893 E-3 O.3I9 E+! 0.105 E-1 0.484 E-2 0.190 E+) 2° 
55.22137 0.6274 E-3 0.262 E+1 0.1079 E-1 0.417 E-2 0.207 E+4 19° 
55.78382 0.9471 E-3 0.212 E+1 0.1110 E-1 0.375 E-2 0.207 E41 WW 
56.26477 0.5453 E~3 0.100 E~1 0.1646 E-1 0.774 E-2 0.890 pi’ 1° 
56.36339 0.1335 E-2 0.166 E+] 0.1144 E-I 0.297 E-2 0.229 Ett 157 
56.96818 0.1752 E-2 0.126 E+1 0.1181 E-1 0.212 E-2 0.253 E+l 37 
$7.61249 0.2825 E-~2 0.910 0.4221 E-1 0.940 E~3 0.376 E41 t- 
58.32389 0.2369 E~2 0.621 0.1266 E-1 -0.550 E-3 —O.11E E42 D2 o 
58.44658 0.1447 E-2 0.827 E-I 0.1449 E-1 0.597 E-2 0.790 3° 
59.16422 0.2387 E-2 0.386 0.1319 E-I —0.244 E-2 0.700 E-1 T 
$9.59098 0.2097 E-2 0.207 0.1360 E-1 0.344 E-2 0.490 5* 
60.30604 0.2109 E-2 0.207 0.1382 E-1 ~0.435 E~2° 0.680 D3 5” 
60.43478 0.2444 E-2 0.386 0.1297 E-1 0.132 E-2 —0.120 E41 7 
61.15057 0.2486 E-2 0.621 0.1248 E-I -0.360 E-3 0.584 E+1 9° 
61.80016 0.2281 E-2 0.910 0.1207 E-1 —0.159 E~2 0.286 E+1 We 
67.4¢122 0.1919 E~2 0.126 E+1 O.tI71 E-1 -0.266 E-2 0.226 E+ D4 13° 
62.48626 0.1507 E-2 0.827 E~I 0.1468 E-I -0.503 E-2° 0.850 3 
62.99797 0.1492 E-2 0.166 E+1 0.3139 E-1 -0.334 E-2 0.218 EFI 1S* 
63,56852 0.1079 E-2 0.212 E+! 0.1108 E-I —0.417 E-2 0.196 E41 17° 
64.12778 0.728! E-3 0.262 E+! 0.1078 E-I -0.448 E-2 0.200 E+} 19° 
64.67886 0.4601 E-3 0.319 E+1 0.105 E-I -0.515 E-2 0.184 E+1 2° 
65.22412 0.2727 E-3 0.38! E+I 0.102 E-I -0.507 E~2 0.192 E+1 23° 
65.76474 0.152 E-3 0.448 E+] 0.100 E-1 —0.567 E-2 0.178 E+} 25° 
66.30195 0.794 E-4 0.522 E+1 0.97 E-2 -0.549 E~-2 0.184 E+] 27° 
66.83663 0.391 E-~4 0.600 E+1 0.94 E-2 ~0.588 E-2 0.174 E+1 29° 
67.36933 0.181 E-4 0.684 E+1 0.93 E-2 -0.560 E-2 0.177 E41 3° 
67,90051 0.795 E-5 0.774 E+1 0.89 E-2 —0.580 E-2 O.I7I Ett 33° 
68.43054 0.328 E-5 0.869 E+1 0.87 E-2 -0.570 E-2 0.165 Et 35° 
68.95972 0.128 E-$ 0.969 E+1 0.86 E-2 -0.530 E-2 0.174 E+1 37° 
118.75034 0.9341 E-3 0.000 0.1592 E-1 -0.441 E-3 0.890 I 





Yable f. (continued) 











Center Quantum number 
frequency, Strength, a, Temperature Width, a, = WentMicatin a). <. 
(GIlz) (kliz/kPa) exponent, a, (Gllz/k?Pa) Lower Upper . 
368.498350 : 0.679 E-4 0.200 E-1 0.156 E-t 1,1 2,3 
424.763120 0.638 E-3J 0.122 E-I 0.147 E-t 2,3 2,3 
487.249371 0.235 E-3 0.122 E-I 0.147 B-} 2,1 3,3 
715.393150 0.996 E-4 0.891 E-1 0.144 EI 3,3 4,5 
773,839732 0.571 E-3 0.798 E-1 0.140 E-1 4,3 4,5 
834.145790 0.180 E-3 0.798 E~I 0.140 E-1 4,3 5,5 
Identification (th,O) 
V% b, b, b, Lower Upper 
. 22.235080 0.105 0.214 E4+I 0.281 E-t 5,2, 3 6, 1.6 
68.052 0.180 E-2 OBIS EVI 0.280 E~I 3,2, 11% 4,1.3 
. 183.310091 0.238 Est 0.653 0.282 E-| 2,2, 0 3,1,3 
32%.225644 0.460 E-1 0.616 Lad 0.220 E-1 9,3, 6 10, 2,9 
© 325.1529t9 0.155 E41 0.182 el 0.290 E-I 4,2, 2. 5,1, 5 
*  380.197372 0.123 E+2 0.102 Ev] 0.285 E-I 3,2, 1 4,14 
, 386.778 0.400 E-2 0.733 Eat 0.160 E-1 $1, 2, 10 10,3, 7 
. 437.34667 0.630 E-1 0.502 Et) 0.150 E-1 6,6, 0 7,5,3 
439. 150812 0.921 0.356 Et) 0.175 E-1 5,5, 9 6, 4,3 
443. 018295 0.191 0.502 E41 0.148 E-I 6,6, | 7,5,2 
*  448.001075 0.107 E+2 O.N37 Eel 0.246 E~1 3,3, 0 4,2,3 
470.888947 0.328 ; 0.357 Ev 0.181 E-1 5,5, 1 6,4, 2° 
474.689127 0.124 E+! 0.234 Etl 0.210 E-3 4,4, 0 5.3.3 
‘ 488491433 0.256 0.281 E14 0.222 E-t 7,1, 7 6,2,4 
, $04.219 0.380 E-1 0.669 E4t 0.127 E-1 71,7, 0 8, 6,3 : 
505.126 0.120 E-1 0.669 Et] 0.130 E-I 7.7 0 8, 6,2 
*  $56.936002 0.526 E43 0.814 0.317 E-I 1,0, 3 1, 1,0 
*  620.700807 0.521 E+! OME 0.216 E-t 4,4, 1 5,3, 2 
658.340 0.460 0.776 Ett 0.328 E-1 10 Iq) 1,1,0 
*  7$2.033227 0.259 E43 0.336 0.302 E-I 2,0, 2 zt 
836.836 0.120 E-1 O.BIf Ett 0.170 E-t 1,2, 9 10, 5, 6 
859.810 0.150 E-I 0.799 Ett 0.270 E-1 2,0, 2(1) 210 
899.380 0.910 E-1 0.784 Ed 0.300 E-4 tl tq) 2,0, 2 
903.280 0.640 E-5 0.835 Ev} 0.280 E~-1 2,2, 11) 3, 1,2 
907.77) 0.179 OSOIE Vt 0.204 E-1 8.3. 5 9,2, 8 
© 916.169 0.890 E+ t O17 EE 0.249 E-t 1,3, 1 4, 2,2 
970.320 0.940 £4) O64 El 0.246 E-I 4,3, 1 5,2,4 
. 987.910 0.145 Ea3 0.180 0.299 E-! Lt, ot 2,0, 2 
: + 1097.368 0.840 E43 0.656 0.335 D3 3,0, 3 3,1,2 








Read 0. 307, E-Jas 0. 307 x 107’. 

“Nonresonant O, spectrum (equation (21)). 

*DI denotes doublet. 

"Rosenkrant’s [1975] first-order solution for the 60-GIiz Aan shape has a shortcoming: the listed values of a,(K = 3°, 


$-) have to be reduced 5% (which is of negligible consequence) to assure that a(O,) 2 0.000 for f > 160 Gilz even when 
N72 O (equation (21). 


“Here (1) denotes first vibrationally excited state. 














Ill. 


SHELLMAN - MORFITT ROOT FINDING ROUTINE 


A. GENERAL OUTLINE 


The Shellman-Morfitt root finding routine was developed by C.H. Shellman and 


. D.G. Morfi 


ionosphere 


ttto locate ELE VILE LP mode constants of wave propagating in an earth- 


Waveguide environment. The theory of the algorithm described in this 


chapter is based on References | and 14. 


The complex root finding routine is used to find the zeros of the modal function. 


The routine 


will locate all simple zeros of an analytic function in a prescibed rectangular 


rezion of the comple\ 2-plane. The principle of the root finding method is based on the 


following theorem [Rer. 15}: 


Theorem (Argument Principle): 

Let ID be a simply connected domain and fiz) be analitic in D except at a 
finite number of poles. Let F be a closed contour in D not passing through 
anv of the zeros or poles ef (7). Then, the accumulated phase change. A @, 


of f(7) around F traversed in a clockwise direction is given by 
Ad = 22(\,— \)). (240) 


where .V, is the number of zeros and NV, is the number of poles (counting 


muluplicities) enclosed by T. 


From the above theorem, if f(z) has only zeros in the region enclosed by T . then 


every constant phase contour associated with a zero crosses PF. 


An analvtic function f(z) can be expressed as 


. where @ is the phase of fiz) given bv 


f(a) =. (Re Y + Unify expiay (241) 

















If {(7) does not have anv pole or branch point, then the constant phase curve. @ = 0.. 
radiating from a zero of f(z} must cross the closed contour, F, enclosing the zero at least 
once. In addition, no other zero of (7) may be on this phase curve. This is illustrated in 
Figure 5. 

A line of constant phase which crosses the contour T may be followed until it leads 
to a zero or until it crosses the contour again. A zero of {(z) can be determined from the 


intersection of the curves 


Im(f) = 0 » (0, = 0 or z) (243) 


and 





Re(f) = 0 2 (0, = = or -) (244) 


tay 
to 


94 


Closed contour [ 





Figure 5.) Constant Phase Line of f(z) [From Ref. 
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B. ROOT SEARCH PROCEDURE 
A search rectangle is set up in the complex z-plane so that the simple zeros of fi7) 
will be located by the routine. The corners of this rectangle are denoted by: 


® ty 


value of the real part of z at the left edge of the rectangle. 
© f,- Value of the real part of 2 at the right edge of the rectangle. 


e 1, - Value of the imaginary part of z at the top edge of the rectangle. 


e 
~ 
7) 
' 


value of the imaginary part of z at the bottom of the rectangle. 


The search rectangle 1s divided into small mesh squares of length 6. Restrictions on 
the size of 6 will be addressed later. 
A new search rectangle is generated from the specified search rectangle by express- 


ing the edge of the new rectangle in terms of mesh units. One mesh uit is of length 6. 


Thus 
ji kes im (4) (248) 
Lp = 5 wet 
! 
Jr = Int i) F (24) 
t- 
Jp = Int (+) (247) 
and 
! 
J_ = Int (<) (248) 


where [nt(x) denotes the integer part of x. In order that the original search rectangle falls 
completely within the new search rectangle, the new rectangle is made one mesh unit 


larger on all sides. That is, 


on Penn th > oa 

= fre 5 L= 16 

J, = LJ _? aoe ee (249) 
2 > 0 

j oni rt Sfp. = ist 

as + | ip ely. * o 
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(tek 2 ae 


= ie 
of ies Lr + | > ip SO oy 
and 
i i ae a (282) 
i ee > tp < 0 


The search rectangle with mesh grid is shown in Figure 6. The pius and minus signs 
at the corners of the mesh squares are the signs of Im(1). The solid curves represent the 
phase line, Im(M= 0, while the dashed curves represent the phase line. Re(f)=0. The 
zeros of {(Z) are located at the crossing of the Im(f)=0 and Re(M=0 curves. A local 
mesh coordinate system is shown in Figure 7, The lower left-hand corner of the mesh 
square is taken as the origin of that mesh square’s local coordinate svstem. The value 
of f(z) evaluated in the local mesh coordinates of the A” mesh square is denoted bv 


Sy, .c,). Thus, the values of f(z) at the corners of mesh square | are 


GOW A ef Oe ETI (253) 
£ Oo, O) = ffl, d + f(Jp —- 6], (254) 
FOU) = fide ONS BIGGEST) (288) 
and 
SPUN = fl, + Dd+jrd]. (256) 


A basic assumption made on f(z) is that Im(f) and Re(f) are linear functions of the 
local mesh coordinates, y, and@, . This means that in the mesh square k. f(z) can be 


expressed as [Ref. 1] 


f Ge ©; ) = U). + hb, Op, + Ce Xk + dj, G), an ‘ (257) 


A) 
‘am 











Figure 6. — Illustration of search rectangle with mesh grid [From Ref. I] 


where 
a = £00), (258) 
yh = f™oOy—-sOjw,o , (259) 
Gq = fo) — £0, 0) (260) 
and 
d& = fPO9)+/7 RH -£%O0) -£O Uo . (261) 


This assumption puts an upper bound on the size of the mesh unit, 6. The mesh 
unit should be sufficiently small to censure that the linear variation of f(z) holds. How- 


ever, 6 should not be excessively smal! to reduce computer run time. 














Figure 7. Local mesh coordinate system [From Ref. 1] 


The zeros of {{z) are searched by examining the imaginary part of {{z) counter- 
clockwise around the rectangular contour. For the example shown in Figure 6, the edges 
of the search rectangle are examined for crossing of the phase curve, Im(f)=9, starting 
with mesh square I, edge I. The values Jin[ f (0, LJ and JmUf (0, 0)J are first com- 
puted. If they are of the same sign, then there will be no Im(f)=0 curve passing through 
edee I. 

Mesh square 2 is investigated next. The values of Jinf f @ (0, IJ and fal fi ® (0, OF) 
are computed. The value of Jin f @ (0, 1)] is identical to Jnl f (0, 0)] . If the signs 
of Im f (0, 1)) and Imp f@(0,0)] are opposite as shown in Figure 6, then the 
Iinn(f}=0 curve enters mesh square 2 along cdge J. The crossing point is given by [Ref. 
I 


yy = 0 (262) 








7 Int f © (0, oY) 
Inf f ©, 1)] — Inf f® 0,99) 


OE (263) 


The next step is to determine if a zero is Within mesh square 2. The values of 


ImC f © (1, 0)] and Dnlf @ (1, 1)] are first computed. Two tests must be made: 


l. 


A test is made to determine if there is one or two Im(f)=0 curves entering and 
leaving the mesh square. If Jm[f  (0,0)] and Jn[ f ® (1, 1)] are of the same sign, 
e.g., minus, and Inf f (0, 1)] and Jnl f (1, 0)] both have opposite signs, c.g., 
plus, then there are two lines of Im(f}=0 entering and leaving the square. This 
condition is depicted in mesh square 14 of Figure 6. Otherwise. there is only one 
line of Im(f }=0 entering and leaving the mesh square. This is illustrated in mesh 
square 2 of Figure 6. 

A test has to be performed to determine if there is at Ieast one Re(f)=0 entering 
and Jeaving the mesh square. If Re[ f (0,0)] , ReLf (0. 1)] , Ref (1. 0)] and 
Re[ f (1, 1)] are all of the same sign, then there is no line of Re(f)}=0 passing 


through the mesh square. Otherwise, there is at least one Re(f)=0 curve entering 
and leaving the mesh square. 


Note that the above tests will fail if there are two or more Im(N=0 or Re(N=0 


crossing the mesh square at the same edge. As shown in Figure 8, two Im(f)=0 crossing 


the same edge of a mesh square results in a zero being missed. In this situation, the mesh 


size has to be reduced in order to resolve the two zcros. 





Perceived Perceived 


Im(f)}=0— \y 





Actual 
Im(f)=0 


wat 








Figue 8. Crossing of two Im(f) = 0 at the same edge of the mesh square 


‘The program (sce subroutine [zerox in Chapter V) also checks for the presence of 
more than one zero on the same phase line. Such a situation may occur as in Figure 9. 
‘Lhe program treats such an occurrence as an error and will proceed to reduce the mesh 
- size. The root-search will start all over again. 


One can infer that the presence of more than one zero on the same phase line is duc 


: to the inability to resolve two phase lines crossing a mesh square at the same edge. These 








zeros actually lie on different phase lines and are ‘legitimate’ zeros. It is the author's 
opinion that these zeros are acceptable and a mesh-size reduction accompanied by a re- 


run of the root search routine unnecessary. 


Perceived 
Re(f)=0 


Re(f)=0 Perceived 
Im(f)=0 





Figure 9. Two zeros found on the same phase line 


‘The line Im(f)=0 is traced from mesh to mesh until it exits the contour rectangle. 
For the example in Figure 6, the curve AB is followed through mesh squares 2, 3, 7, 1, 


10, 14 and 13. In mesh square 11, both Im(f)=0 and Re(f=0 pass through the mesh 
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square. The point of imtersection can be determined from Equation 287. The math- 
ematical details are presented in section C. 

The contour AB is followed out of mesh square I] into mesh squares [0 and I4. 
In mesh square 14, there are two Im(f=0 curves and one Re(f)=0 curve. A crossing 
point is chosen to be a zero of f(z) if it hes on the Im(fM=0 curve currently bemg fol- 
lowed. Therefore, when the contour AB is being traced. the zero in mesh square 14 will 
be ignored. Hlowever, this zero will be located when contour CD ts being traced. 

When contour AB exits the search rectangle at point B. the value of the leading 
edge (counterclockwise) of the mesh square is noted. In this example. the value of the 
leading edge is f (0, 1). Identifving the exit points would prevent subsequent re-enter- 
ing of the search rectangle via these points. 

After the contour AB exits at B, the search for zeros of f(z) 1s continued by looking 
for more crossings of Im(f)=Q with the search rectangle. Since mesh square 2 was the 
Just mesh square that was checked for such crossings, the next mesh square to be ex- 
amined will be mesh square 3. 

The next crossing of Im()=0 with the search rectangle is at mesh square 12. This 
contour CD is again followed through to mesh 21 where the exit point D is noted. The 
zero along CD in mesh square 14 is considered as an acceptable zero of f(z) since it hes 
on the contour CD which is the current contour being traced. The rectangular contour 
search will continue as described above until the search gets to mesh square 1, edge 4 


where it stops. 


C. DETERMINATION OF THE COORDINATES OF ROOTS 
The root of f(z) 1s the point of intersection of Im(M=0 and Re(A =. This can be 
determined from Equation 257, 


By equating Im(f) of Equation 257 to zero, the following equations are obtained: 


r Ima.) + Imth,eo, (264) 
a Lin(e,) + Lin(d berg ara 


and 


7. Finlay) + Irate Vy, (263) 
os a Im(b,) + Inka 7): oe 


Similarly, by equating Re(f) to zero, 


6] 


Rela.) + Rethy eo; re 
Le = Rete.) + Rel deo; oa 


and 


Rel a.) + Ref Ch Zp 


O, = Reth,) + Redz, 


(267) 


Note that Equations 264 and 268 are equivalent; so are Equations 266 and 267, 
By equating the right hand side of Equation 264 and that of Equation 266, the following 


quadratic equation in w, results: 
Pyeoy” + Qo, + We = 0 ’ (268} 
where 


r, = Re (h,.) Ind, =< Tn( by) Re (ci) ; (269) 


Q, = Re (ay) fin(d.) + Re (h,) Ini(e;) 


(270) 
— danila) Re la.) — Inih,) Re (ez) 
and 
W, = Re (a,) lnile,) — Inila,) Re ley). (271) 
Similarly, from Equation 265 and 267, the quadratic equation in 7, is obtained: 
5 es + OQrv7, + Hy = 0, (272) 
Where 
Py = Re(cy) Ind) — Inic,) Rela), (273) 
Q, = Relay) Ind.) + Rele;) lin(b,) 
(274) 
— Ima,) Re(d,) ~ Tin(c,) Re(h,) 
and 
Wy = Rela, inth,) — Tnita,) Reth,) . (278) 


According to the program source code, if 
IPL <1 Ps (27% 


then the root-finder routine uses Equations 264 and 268 to solve for y, andw, . If the 
condition m Equation 276 is not satisfied. then 7, and@, will be determined from 
Equations 2645 and 272. 

Hlowever, it should be noted that if Equation 276 is true and if | ?, | approaches 
zero, one of the solutions of Equation 268 will become infinite. By the same token. if 
Equation 276 is not true and if | ?,| approaches zero, one of the solutions of Equation 
272 will become infinite. The program (see subroutine quad in Chapter V) prevents such 
occurrences by incorporating a test such that the program returns only one solution 
when | ?,} or | P| approaches zero, 1e.. Equation 268 or 272 reduces to a linear 
equauon. 

On the other hand..unless both | 7)| and | P,] become very small, the above 
problem could be avoided if 7, and «, are determined from Equations 265 and 272 when 
Equation 276 is true: otherwise Equations 264 and 268 are used. This modification is 
reconumended. 

A solution of Equations 265 and 272 or 264 and 268 will be valid if it Hes within the 


current mesh square, Le. 


(273) 


where zc, = (7, @.) 1s the solution . 
In addition, the solution must also lie on the line Im(f}=0 currently being traced. 
The test for this condition is derived below. 


From Equation 257, Im(=0 and Re(f)=9 are equations of an equilateral hvper- 


bola with vertical and horizontal asymptotes given by 


In (bh) 


te Im (dd) 





and 


Tine) 
Oe = Sens (280) 
Let (77, @;) be the point at which Im()=0 enter the mesh square. If 
Wee Jel ea Sod (281) 
and 
(Of ~ Oc) (MW, —- @) > a. (282) 


then the solution of f(z) lies on the current line. 
Figure 10 illustrates the case where 2, is a proper solution and Tigure 1] shows the 


cuse Where 2, is not a proper solution of [(z)=0. 
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(Xe ’ W,) 
(Xe ’ W) 


f(0, 0) f(1, 0) 


Figure 10. z, is a proper solution of {(z)= 0 [After Ref. 14] 





(X., w,) 





ee EEE SEEGER 


ligue IJ. c, is not a proper solution of {(7) = 0 {After Ref. 14] 
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D. SOLUTION OF THE QUADRATIC EQUATION 
Procedures for the soludion ef the quadratic equavon are required for solving 
Lyuauony 268.and 272. 


Consider the quadratic eqguauen of the ferm 


ax + 2hbxy to = 0, ( 


tu 
nm 
ae 


The solution is 











-f[-1: ee =| (284) 
Bt i 
Se ease 
Where 
b= 2c (288) 
Hence. 
. w= of-14+ Tee] (286) 
and 
yy = a [-t-\t+e] , (287) 


Ife << [. then 








Thus, 


1 | ! 
+(4-1)(4- 
ac 


A —ae = F 
Factoring out (>) . Equation 289 becomes 
i 


(290) 


(291) 


In general, 





and 


ble SEE ae ee ; 


205 
h 


Once x, is found, x, can be found from Equations 286 and 287. By adding Equations 


286 and 287, x, is given by 
(296) 


E. NEWTON RAPHSON ITERATION 
The zeros of {(7) found by the Shellman and Morfitt routine are only approximate 
solutions. The Newton Ruphson iteration is used to improve the accuracy of the zero 


locations. The new root, 2... 1s related to the previous root, 2. bv 


(297) 


(298) 


Iterauion of Equation 297 is continued until Az, is smaller than a pre-assigned tolerance. 
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IV. SOLUTION OF STOKES’ DIFFERENTIAL EQUATION 


A. STOKES’ EQUATION 

When the modified refractivity profile is modelled by a series of piecewise linear 
segments in the mathemaucal formulation of wave propagation in a mululaver tropos- 
pheric waveguide environment, the height gain function in each layer can be obtained 


from solution of Stokes’ Equation. 


a2 
A 


Se oe GP Ap la =O : : 5 (299) 


Cd 
Solutions of Stokes’ Equation are commonly given in terms of Airy functions or in terms 
of modified Hankel functions of order one-third. Several possible solutions exist and the 


choice of solution is dependent on numerical and physical considerations. Since Ji{q) 


is proportional to the attenuation rate of the electromagnetic signal, physical consider- 


auons dictate that Jui(g) =O. In ths program. the height gain function is expressed 
as a linear combination of A, (gy) and &,(qy) which are linearly independent for 
Imig) 2% 0. VYhe choice of A, (g), Ar (g) or equivalently {4i( — ge U8). A= gy df as 
to assure that at least one of the two functions do not become very small numerically in 
the upper half of the complex g-plane. For the top layer, i.¢.. i= 1. the boundary condi- 
tun at 2 > oo requires that the height gain function represents an outgoing wave. 
Hence in the topmost laver. A, (y,) or equivalently -f(— g,e~?7+) is chosen to saust¥ the 


radiation condition. 


B. EVALUATION OF AIRY FUNCTIONS 
I. General Outline 

The algorithm for evaluation of complex Airy functions is adopted from Z. 
Schulten and D.G.M. Anderson [Ref. 16]. Evaluation of the complex function over the 
entire 7-plane was achieved with a Taylor series and generalised Gaussian quadrature 
method. A Tavlor series is used only for a smal] region near the origin within an ellipse 
with foci (-1.58, 3.95) and (0.28, —2.11), major axis of 9.0 and munor axis of 6.4. 
Outside the ellipse, the Gaussian quadrature approximation is used. The elliptical 


boundary was chosen to satisfy restricuons imposed by computauonal methods. 





Analysis of the Airy funcuons need only be done in the upper half complex 


plane as those in the lower half plane can be computed by the following conjugate 


relations: 





dAi(z) dfi (2) 
—_—_——- = ——_., (309) 
dz de 


Therefore, an algorithm can be developed to compute Antz) for O < argz < 22/3 and 
by conjugacy, An?) for -27/3 < args < Ocan be computed. 
For 22/3 < arg: < 42/3, the following connection formula of Airy functions 


can be used: 


dite see ae ae (302) 


2. Evaluation by Taylor Series Expansion 
Por small] values of |7 |, the Airy functions are evaluated by the series expansion 


[Ref. 1 and 4}. 


and 


AM(2) = Of’ (2) — Coe’ ts). (304) 


where 
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(3060) 


= - r (k++) _3k+1 
pk gree es ees 307 
& (>) : 2) Gk+Ip (307) 


g'(:)= )3 ——— =— .. (308) 


; ayer} 2) —-) 
Cc, = Ee r(+)] (309) 


and 
one eres) (310) 


The number of terms required for convergence of the series is determined by a tolerance 


factor. e. such that 


n~] 
v 
[a] om e Pa (311) 
i=0 
where 
n 
Ai) = Sa G12) 








Similiar test for convergence is unplemented for /(2). In the program, ¢ is taken to be 


the machine epulen, Le. 


Evaluation by Generalised Gaussian Quadrature Method 
given by [Ref. 1¢] 


2 


Integral expressions for A(z) 1 





a 
: J “Vig —2 plxidx 
Ai(z) = oe _ 
a Len 
™ mu (314) 
ye 
for jel Oy. dates <= 
J 
where 
’ 2 sak heoos 
s= S27 (ois) 
2? 
and 
ee re Meee ee re Ay ‘ 
AX) = x 2 3 x ¢ Ai( 7 ) : (316) 


The function, p(x). is non-negative and exponentally decreasing. The integral in 
Equation 344 can be evaluated by the generalised Gaussian quadratic approximation, 


Le, 


”n n 
pix) dv Wi ss 
ee os (317) 
ép » +N ra aie Xj 
‘ =i 


Where w. is the quadrature weight and x, is the abscissa. Hence 


a 
. baie ee See NS 446 
AME) SN Se El 318 
(>) 5 a ree (318) 








Z oe 1 
aie | es roses Po ee Moy Wi ¥ Wy Bee 
Ai(c) ~ ze Wea e > =. ———_ -: ——_— . (319) 
2 6 TP Zhi tx) . e 2 
— Sea a4 (. +.x;) 
i=1 i=) 


The quadrature weights, w, . and abscissac, x, had been calculated by Z. Schulten et al. 
[Ref 16] and are given in Tables I] to TV. Lhese tables had been used for computation 


of single precision values up to seven significant digits. 
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Table 2. 


Table 3. 


NUMBER OF TERMS REQUIRED IN QUADRATURE FORMULA 


[AFTER REF. 16] 


z arg z Ai(z) 
2a 
jz) > 11 <s 2 ferms 
2a 
§5<[z|< 1] i 4 terms 
vis 
2.5<|z|]<5 = 4 terms 





4-TERM GENERALISED GAUSSIAN INTEGRATION FOR AIRY 


FUNCTIONS [FROM REF. 16] 


x4 


3.9198329554455091 
1.6915619004823504 
5.0275532467263018(—01) 
1.9247060562015692( —02) 





We 


4.7763903057577263(—05S) 
4.99 14306432910959(—03) 
8.6169846993840312(—02) 
9.0879095845981102(—01) . 








Table 4. 2-TERM GENERALISED GAUSSIAN INTEGRATION FOR AIRY 
FUNCTIONS [FROM REF. 16] 


Xj W; 


1.0592469382112378 3.1927194042263958(—02) 
3.6800601866153044(—02) 9.6807280595773604(—01) 


C. EXTENDED COMPLEX ARITHMETIC 

The extended complex arithmetic was introduced by Baumgartner in NWYVG to 
handle complex numbers of large magnitude to avoid overflow problem. 

It was discussed in Section D of Chapter IH that for large values of |g |, the mag- 
nitude of 4,(g) and k,(g) may become exponentially large or exponentially small. 

Numerical evaluation of the modal function of Equation [66 can easily vield com- 
plex numbers with magnitudes as large as 10°! or as small as 10°! . Numbers of these 
magnitudes are outside the numerical limits of most computers. 

To overcome this problem, each complex function, &,(g) (or k,(q)), is represented 
by 


ce Ee 5 (320) 


Where k, is the complex amplitude and @ is the real exponent. In the program, 4, is 


“normalized” such that 


eT} max Re (hy) |. Lon) 1 (321) 








and ¢ has integer values. By evaluating the complex amplitude and real exponent sepa- 


rately, overflow or underflow problems are avoided. Functions expressed in the form 
as in Equation 320 are known as extended complex numbers. Numerical manipulation 
such as addition and multiplication involving extended complex numbers are performed 
in the following manner. If z, = Z, e¢: and z, =2,e*2 are two extended complex num- 


ber, then their product is 


zZ = 2; 2, 
ate, (322) 
where 
ZF = 2, 2, (323) 
and 
¢=O,+62 . (324) 
The sum of two extended complex numbers is given by ; 
where 
Batic) 5 6 > by | (326) 
and 
¢$=6 3¢2 & (327) 
or 
P= 2+ 2,e%-%) 3 g < (328) 
and 
@ = $ 3 O< 2 - (329) 
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V. DESCRIPTION OF MLAYER PROGRAM ELEMENTS 


A short description of the functions of each program clement in MLAYER is pro- 
vided in this Chapter. The program had been tested and ran successfully on the Micro- 
solt FORTRAN version 5.0 compilers and Microsoft C compilers version 5.0. The 
maximum number of lavers and modes are 35 and 150 respectively. These limits are the 


upper bounds for operation on personal computers with 640 kilobytes of main memory. 


A. MAIN 
i. Description 

Main is the controlling program clement for MLAYER. It calls appropriate 
subroutines to calculate modal eigenvalues, mode sum, path loss and radio horizon dis- 
tance. Main also calls subroutine wvgstdin to read in data from input file, filcin. Three 
output files are opened, namely, filein.out, filein.eig and filein.plt. 

 Filein.out provides a documented output of input parameters and a list of mo- 
dal eigenvalues found in each search rectangle. Filein.out also contains a list of eigen- 
values in increasing order of real parts with the corresponding values in the @-plane and 
respective attenuation of cach mode in decibels per kilometer. 

Filein.cig is an exact replica of the input file except that the parameter, myfile, 
is set to one and it contains a list of eigenvalues found. Once filein.eig is created, it can 
be used as the input file to re-run MLAYER for various transmitter,receiver range and 
height configurations. 

Filein.plt contains cight columns of data, namely range. transmitter height, re- 
ceiver height. coherent mode sum. incoherent mode sum, coherent path loss, incoherent 
path loss and radio horizon distance. The data are arranged in this manner to facilitate 
graph plotting. 

Samples of filein, filein.out, filein.eig and filcin.plt are given in Appendix B. 

In addition, main also computes the following: ; 
e The derivative of the modified refractivity with respect to z for cach layer. 
e The derivative of the atmospheric absorption with respect to z for each layer. 
¢ The absorptivity and its derivative with respect to z for each layer. 
e The derivative of q with respect to z for each layer. 


@ The modal attenuation rate for each mode. 


78 





e The coefficients required for computation of surface roughiness, 


¢ The eigenvalues in the @-plans. The eigenvalues in the @-plane are related to those 
in the q,,-plane by 


eer ae in-F Es p° | (330) 


The conversion to @-plane is done solely for comparison of results with other models 


where root-search was performed in the @-plane. 


2. Inputs 


Input data are read in by subroutine wvgstdin in the following order: 


¢  filein 


e = mifile 

e = fgmzin 
© delfq 

© nfreq 

® mpol 

e aloss 

® seatmp 
e = seaslt 

e iflgab 

®  airtrap 
e rh 

e wepm3 
e rmsbht 
e ztinit 

e delzt 

e nzt 


: Pathname of input file. 


: If mfile=0, read input data and calculate eigenvalues. 
If mfile= 1, read input data and eigenvalues. 


: Initial frev’’cncy in megahertz. 
: Frequency increment in megahertz. 
: Number of frequencies for which modes will be found. 


: If mpol=0, wave is horizontally polarized. 
If mpol= 1, wave is vertically polarized. 


: Maximum rate of attenuation in decibels per kilometer of modes 
that will be found. 


: Scawater temperature in degrees Celsius. - 
: Seawater salinity in grams of salt per kilogram of seawater. 


: If iflgab=0, calculate the height invariant atmospheric absorption. 
If iflgab= 1, read in array containing values of atmospheric 
absorption at laver boundaries. 


: Air temperature in degrees Celsius at reference height, Zref. 
Airunp is used for calculating the height invariant atmospheric 
absorption when iflgab= 0. 


: Relative humidity in percent at reference height. Used for calculating 
the height invariant atmospheric absorption when iflgab=0. 


: Water concentration in grams per cubic meter at reference height. 
Used for calculating the height invariant atmospheric absorption 
when iflgab= 0. 


: Rms bump height in meters of the sea surface. 
: Initial height of transmitter in meters. 
: Height increment of transmitter in meters. 


: Number of transmitter heights at which field strength is to be 
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computed. 


zrinit : Initial height of receiver in meters. 
delzr : Height increment of receiver in meters. 
nzr : Number of receiver heights at which field strength is to be computed. i 
xinit : Initial range separation of the transmitter and receiver in kilometers. 
delx : Increment in range separation in kilometers. ‘ 
nx : Number of range separation values at which field strength is to be 
computed. 
zref + Reference height at which air temperature, relative humidity and 
Water concentration are specified. 
nzlavr : Number of linear segments used to model the modified refractivity 
profile. 
Zi(j) : Array containing values of heights in meters at layer boundaries. 
zimi(j) : Array containing values of modified refractivity at layer boundaries. 
zigab(j) — : Array containing values of atmospheric absorption in decibels 
per kilometer at laver boundaries. 
3. Outputs ° 
The outputs to filein.out (logical file 16) are: 
nrmode =: Number of modes found. : 
qeigen(j) : Array containing values of modal eigenvalues found in the complex 
M.a-plane. 
theta : Complex eigenangle referenced to ground level for cach mode. 
atnu : Rate of attenuation in decibels per kilometer for each mode. 


The outputs to filcin.eig (logical file 17) are : 


nrmode =: Number of modes found. 
qeigen(j) : Array containing values of modal eigenvalues found in the complex 
ia-plane. 

The outputs to filein.plt (logical file 15) 
ecms : Coherent mode sum field strength relative to free space in decibels. ; 
eims : Incoherent mode sum field strength relative to free space in decibels. 
ecpl : Coherent mode sum rath loss in decibels. . 
eipl : Incoherent mode sum path loss in decibels, 
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4. Calling subroutine 


® None 


3. Subroutines called 


© wvestdin 


e chkmod 

® seah2o ; commented out 
e fndmod 

®  casin 

® «ao2h2o ; disabled 


® modsum 
® dhoriz 


¢ function ibstrip 


The call to subroutine seah20 was commented out in revision 3.0 of the pro- 
gram. To use the subroutine, remove the comment statement. Subroutine ao2h2o was 
disabled by setting iflgab to one. To enable ao2h2o, set iflgab to zero. 


6. Common block areas 


¢ com |! 
® com 2 
® com 3 
e pap | 

e datum 


B. SUBROUTINE ABCOEF 
I, Description 
Subroutine abcoef calculates the coefficients, 4, and B, of the height gain 
function given by Equations 25 and 28. Computations of the coeflicients are done in 
extended complex arithmetic. The normalizing integral of Equation 77 is also evaluated. 


The coefficients are then normalized by the following relation: 
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By noting that 


5 


=f py 





it can be seen that the denominator of Equation 331 is identical to that of Equation 77. 
In calculating the coefficients, the decision to integrate up or down is deter- 


mined by the rate of attenuation given by 


att = Re[pp] x20 loge x 10’ dBikm 


If the attenuation rate is greater than 0.1 dB km, upward integration is performed; oth- 
erwise downward integration 1s carried out. 
At the same ume, the subroutine also checks for evanescence using Equations 
153 through 156 for upward integration and Equations 162 through 165 for downward 
integration. 
2. Calling statement 
The subroutine is called by the statement 
Call abcoef (zero, m). 
3. Inputs 
The input variables are 
zero : An eigenvalue in the q,,-plane . 


m : Mode number. 


4. Outputs 


The output variables are 


acoefa(i.m): Two dimensional array element of complex amplitude of A in the 
@ Javer and or" mode. 


acoefe(im) : 1wo dimensional array element of real exponent of A in the 
” Javer and nr mode. 





¢ = beoelugimy: Two dimensional array element of complex amplitude of B in the 
® laver and we mode. 
¢ beoeletim): Two dimensional array element of rea! exponent of Bin the 


i laver and ve mode. 


5. Calling program element 


¢ modsum 


6. Subroutines called 
® xedaj 


® xcadd 


7. Common block areas 


* com | 
* com 2 
¢ pap l 
® pap 2 


C. SUBROUTINE ADDX 
1. Description 


This subroutine adds two extended real numoers according to the following 


equation: 
reit¢ 2 (334) 
where 
z= rz exp(r)  . (335) 
Z) = 2g EXP (zy) (336) 
and 
Tp = 2a, eEXP(Za,) . (337) 


2. Calling statement 
The subroutine is called by the statement 


Call addx (za. ze. zla. ze. z2a. 72¢). 





3. Inputs 


The input variables are 


@ zla : Real amplitude of the extended real number, 2, . 
© zic : Real exponent of the extended real number, <.. : 
® Z24 : Real amplitude of the extended real number, -, . 
@ 22¢ : Real exponent of the extended real number, 7, . : 


4. Outputs 
The output variables are 
© 7a : Real amplitude of the extended real number, z. 


* 2¢ > Real exponent of the extended real nuiuber, z. 


§. Calling program element 


® modsum 


6. Subroutine called 


® —normre 


7. Common block area i 


® none 


D. SUBROUTINE AO2H20 
1. Description 
This subroutine calculates the atmospheric absorption coefficient in decibels 
per kilometer due to the quantum mechanical resonances of oxygen molecules and water 
vapour. Computed values are good for frequencies in the range 1 Ghz to 1000 GHz. 
2. Calling statement 
The subroutine is called by the statement 


Call ao2h2o (atmabs, fhz, tc, rh, zm, Wgpm3). 


3. Inputs 
The input variables are 
© fhz : Prequency in Hertz. F 
e tc : Temperature in degrees Celcius. 
¢ orh : Relative humidity in percentage. , 


¢ zm : Height above sea level in meters. 





¢ wepm : Water concentration in grams per cubic meter at reference height. 


4. Output 
Phe output variable is 


atmabs > Atmasphetric absorption coefficient in decibels per kilometer. 


5. Calling program element 


main 


6. Subroutine called 


none 


7. Common block area 


none 


E. SUBROUTINE CHAKMOD 
1. Description 
This subroutine checks for modu] eigenvalues with zero value and discard them. 
A zero-value eigenvalue may sometimes lead to a divide-by-zero error in the mode sum 
calculation. 
2. Calling statement 
The subroutine 1s called Sv the statement 
Call chkmod (qeigen. nrmode). 
3. Inputs 
The input variables are 


e geigen(j) + Array containing location of modal eigenvalues in the complex 
Qaplane. 


® nrmode : Number of modes. 


4. Outputs 
The output variables are 


qeigen(j) : Array containing location of non-zero eigenvalues in complex 
g,,-plane. 


nrmode : Number of modes. 





5, Calling program element 


® main 


6. Subroutine called 


6 
* none 
7. Common block area : 
e none 
F. SUBROUTINE DETR 
1. Description 
This subroutine evaluates the modal function determinant of Equation 47 in 
extended complex arithmetic. The algorithm is discussed in Section Fo of Chapter H. 
2. Calling statement 
The subroutine is called by the statement 
Call detr (deta, dete. nzlavr. d. de. dp, dple, dmi. dmite, dp2. dp2e, din2. dine). 
3, Inputs 
The input variables are 
© onzijayr > Number of linear segments used to model the modified refractivity : 
profile. 
ed : Array containing complex amplitudes of elements along the niain 
diagonal of the modal function determinant. : 
ee > Arrey containing real exponents of elements along the main 
diagonal of the modal function determinant. 
« dp! : Array containing complex amputudes of elements along the diagonal 
one above the main diagonal of the modal function determunant. 
e dple > Array containing real exponents of elements along the diagonal one 
above the main diagonal of the modal function determinant. 
¢ dml : Array containing complex amplitudes of elements along the diagonal 
one below the main diagonal of the modal function deternunant. 
e dmile : Array containing real exponents of elements along the diagonal one 
below the main diagonal of the modal function determinant. 
e dp? : Array containing complex amplitudes of elements along the diagonal 
two above the main diagonal of the modal function determinant. 
e 6dp2e : Array containing real exponents of elements along the diagonal two 
above the main diagonal of the modal function determinant. : 
r ¢ dme2 : Array containing complex amplitudes of elements along the diagonal 
two below the main diagonal of the modal function determinant. 
e dm2e : Array containing real exponents of elements along the diagonal two 
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below the main diagonal of the snodal funcuon determinant, 


4. Outputs 


The output variables are 


e deta : Complex amplitude of the modul function. 
e dete : Real exponent of the modal function. 
5. Calling program element 
© fetvlx 
6. Subroutine called 
@ xcadd 
7, Common block avea 
* none 


G. SUBROUTINE DHORIZ 
1. Description 
The subroutine calculates the radio horizon distance between the transmitter 
and receiver. The radio horizon distance in kilometers is given by 


d= 2. /=2 [ig + yep]. (338) 


\ A 3 oie 


Where a is the radius of the earth, and =, and 2; are the receiver and transmitter heights 
respectively. 
2. Calling statement 
The subroutine is called bv the statement 
Call dhoriz (dhz, zrevr, zxmtr). 
3. Inputs 


The input variables are 


®  ozrevr : Receiver height in meters. 
® zxmtr : Transmitter height in meters. 
4. Output 


The output variable is 


e odhyv : Radio horizon distance i kilometers. 














5. Calling program element 


e main 


6. Subroutine called 


® none 


7. Commom block area 


® none 


H. SUBROUTINE DXDETR 
1. Description 
This subroutine evaluates the modal function determinant and its derivative in 
extended complex arithmetic. Values of the modal function determinant and its deriva- 
live are required by the Newton-Raphson root-finding routine. The methods for evali- 
ating the modal function determinant and its derivative are discussed in Section E of 
Chapter IT. 
2. Calling statement 
The subroutine 1s called bv the statement 
Call dxdetr (deta, dete, dxdeta, dxdete, nzluvr, d, de. dp}. dple, dmi, dmle. dp2. dp2e. 
dm2. din2e. dxd, dxde. dxdpl. dxdple. dxdml. dxdmle. dxdp2. dxdp2e, dxdm2, dxdin2e). 


3. Inputs 
The input variables are . 
e onelavr : Number of linear segments used to model the modified refractivity 
profile. 
e od : Array containing complex amphtudes of elements along the main 


diagonal of the modal function deteniinani. 


e de > Array containing real exponents of elements along the main 
diagonal of the modal function determinant. 


e 6dpl : Array containing complex amplitudes of elements along the diagonal 
one above the main diagonal of the modal function determinant. 


e dple : Arrav containing real exponents of elements along the diagonal one 
above the main diagonal of the modal function determinant. 


e¢ dml : Array containing complex amplitudes of elements along the diagonal 
one below the main diagonal of the modal function determinant. 


e dmie : Array containing real exponents of elements along the diagonal one 
below the main diagonal of the modal function determinant. 


e dp2 > Array containing complex amplitudes of elements along the diagonal 
two above the main diagonal of the modal function determinant. 








dpe > Array contunming real exponents of elements along the durconal two 
wbove the main diagonal of tae modal function determinant, 


2 Array contaunung compen amplitudes of elements along the diagenal 


wo below the niin diagonal of the modal function determunint. 


dme > Array contaming real exponents cf elements along the diagonal two 
below the main diagonal of the modal function determinant. 


dxd > Array contanung comples amplitudes of the derivauves of the 
elements along the main diagonal of the modal function determinant. 


dxde DArray containing real exponents of the derivatives of the elements 
along the main diagonal of the modal funcuon determinant. 


dvdp} > Array containing complex amplitudes of the derivatives of the 
elements along the diagonal one above die main diagonal of the 
modal function determinant. 


dvdple > Array contuming real exponents of the derivatives of the 
elements wiong the diagonal one above the main diegonal af the 
modal function deternunant. 


dxdmil Array contaming complex amplitudes of the derivauves of the 
elements along the diagonal one below the mii diagonal of the 
modal function determinant. 


dxdmile : Array contaiming real expenents of the derivatives of the elements 
along the diagonal one below the main diagonal of the model 
funcuon determinant. 


dadp2 > Array containing complex amplitudes of the derivetives of the 
viements along the diagonal two above the main diagonal of the 
modal function determinant. 

dvdp2e > Array containing real exponents of the derivatives of the elements 
along the diagonal two above the main diagonal of the modal 
function deternunant. 

dxgia2 > Array contaiming compiexs amplitudes of the derivatives of tho 
elements along the diagonal two below the main diagonal of the 
modal function determinant, 

dxdm2e —: Array containing real exponents of the derivatives of the elements 
along the diagonal two below the main diagonal of the modul 
function deternunant. 


4. Outputs 


The output variables are 


deta : Complex amplitude of the modal function. 

dete : Real exponent of the modal function. 

dxdeta : Complex amplitude of the derivative of the modal function. 
dadete > Real exponent of the derivative of the modal function. 
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3. Calling program element 
e  fdfdtx 


6. Subroutine called 


® xcadd 


7. Commom block area 


e none 


I. SUBROUTINE FCTVLX 
1. Description 

This subroutine calculates the values of the elements in the matrix of the modal 
function of Equation 47 and invokes subroutine detr to evaluate the metal function. The 
subroutine shifts the imaginary part of ¢,, by qshift so that the real axis will not fall on 
the mesh line of the search rectangle set up by the Shellman-Morfitt routine. This allows 
the zeros on the real axis to be located. It is known that zeros near or on the negative 
real axis contribute significantly to the fields within a duct. 

However, the shift implics that the Shellman-Morfitt routine will search for 


zeros of ||A (z,,)|| instead of |A(q,,)i] where 


hi =. —Jashifi . (339) 


To compensate for the error, the shift back to g,, is done in subroutine fndmod for all 


zeros found. 
2. Calling statement 
The subroutine is called by the statement 
Call fetvlx (ql lin, deta, dete). 
3. Input | 
The input variable is 


e qllin : Location in the complex q,,-plane where the modal function is to 
be evaluated. 


4. Outputs 
The output variable are 
e deta : Complex amplitude of the modal function. 


e dete : Real exponent of the modal function. 
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5. Calling program element 


@ = fzerox 


6. Subroutines called 


© xedui 
e = xcadd 
e § detr 

¢ surf 


7. Common block areas 


* coml 
* com2 
®* com3 
* cons 


J. SUBROUTINE FDFDTX 
1. Description 
This subroutine calculates the valucs of the elements in the matrices of 
ACg,,)i and a |A(g,,)] and then invokes subroutine dxdetr to evaluate the modal 
function and its derivative. The subroutine also shifts the imaginary part of ¢,, by qshift 
as in subroutine fetvlx. It is similarly compensated by a “shift-back” in subroutine 
fndmod. 
2. Calling statement 
The subroutine is called by the statement 
Call fdfdtx (ql lin, deta, dete, dxdeta, dxdete). 
3. Input 
The input variable is 


e = qllin : Location in the complex q,,-plane where the modal function and its 
derivative is to be evaluated. 


4. Outputs 
The output variables are 
e deta : Complex amplitude of the modal function. 
e dete : Real exponent of the modal function. 
e = dxdeta : Complex amplitude of the derivative of the modal function. . 
e = dxdete : Real exponent of the derivative of the modal function. 
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5. Calling program element 


® {rerox 


6. Subroutines called 


® xcdai 
® xcadd 
e  dxdetr 
e surf 


7. Common block areas 


® coml 
®* come 
* com 
* com 


Kk. SUBROUTINE FINDFX 
1. Description 
The incorporation of the effects cf surface roughness into the mathematical 
model results in discontinuity of the moda! function along the imaginary q,,-axis . The 
Shellman-Morfitt routine will fail if the edge of the search rectangle lies on the imaginary 
aXis. 
Subroutine findfx will determine if anv edge of the search rectangle is on the 


imaginary axis. Hf so, a small effser of 1 x 10° is introduced. If the left edge is on the 


imaginary axis, then the modal function is evaluated at Re(g,,) = + effser. Ifthe right 
edge is on the imaginary axis, the modal function is evaluated at Re(g,,) = — offser. 


2. Calling statement 
The subroutine is called by the statement 
Call findfx (jr, ji, f fe, tleft. tright). 


3. Inputs 
The input variables are 
° jr : Real part of the mesh square coordinates in mesh units. 
e ji : Imaginary part of the mesh squares coordinates in mesh units. 
e  tleft : Real part of g,, at the left edge of the search rectangle. 
e  tright : Real part of g,, at the right edge of the search rectangle. 
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4. Outputs 
The output variables are 
ef : Complex amplitude of the modal function, 


© fe : Real exponent of the modal function. 


S. Calling program element 
e §${Zerox 
6. Subroutine called 


e fetvix 


7. Common Ddlock area 


® {mecom 


L. SUBROUTINE FNDMOD 
1. Description 
This subroutine sets up the rectangular region in the complex gy. ,-plane for the 
Shellman-Morfitt root finding routine to locate the zeros of the modal function. 


The initial search rectangie is set up with 


the left edge of rectangle, theft = 0.0, 


e the sight edge of rectangle. tight = test + tstep, where tstep is the length of the 
rectangle. 


e the top of the rectangle. ttop given by [Ref. 4] 


_ 2x10 loss py (423 a 
ttop = k20 loge re ( oy ) | : (340) 


Where Aloss is the maximum attenuation rate in decibels per kilometer of those 
modes to be found. 


e =the bottom of the rectangle, thot = — tol, where tol is the tolerance to which zeros 
are to be located. 

In order that the zeros on the real axis can be located. a small offset of qshift 

is introduced in subroutines fdfdtx and fetvlx so as to avoid having a mesh line on the 


real axis. This offset is corrected in the subroutine fndmod by the statement 


. 


zeros(k) = zeros(k) — jqshift . (3-41) 


After the initial rectangle has been searched for zeros of the modal function, a 
new search rectangle is set up to the left of the initial rectangle. If no zero is found in 


three consecutive search rectangles, the search for roots to the left of the initial rectangle 
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is stopped. Fndmod next forms search rectangles to the right of the initial rectangles. 
If no zero is found in three consecutive search rectangles, the search for roots to thie 
right of the initial rectangle is stopped. 

Since subroutine {zerox extends the search rectangle by one mesh unit on all 
sides, fndmod also checks for and eliminates all zeros outside the current search 
rectangle. 

After all the zeros have been found, fndmod sorts them in order of increasing 
real parts and stores them in the array called zeros. 

There are two versions of fndmod. namely, fndmod20.for and fndmodS$0.for. 
Fndmod$&0 sets up a smaller mesh and takes longer time to find the modes. PndmodS0 
should be used onlv if fndmod20 fails to find all modes. 

2. Calling subroutine 
The subroutine is called bv the statement 
Call fndmod (qeigen, nrmode, dmdz, detadx, zim. nzlavr, aloss, waveno). 
3. Inputs 


The input variables are 


ns d\f 
¢ dmdz : Array containing values of 
di 
ad oe dy 
e = detadz : Array containing values of —— 
da 
* zim » Array containing values of M(z) at laver boundaries. 
@ onziavr : Number of linear segments used to model the modified refractivity 
profile. 
® aloss : Maximum attenuation rate (in decibels per kilometer) of modes to 
be found. 
® waveno : Wave number. 
4. Outputs 
The output variables are 
© qeigen : Array containing the zeros of the modal function. 


nrmode : Actual number of modes found. 
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Calling program element 


main 


6. Subroutine called 


{7erox 


7. Common block areas 
coni4 


coma 


SUBROUTINE FZEROX 

Description 

Subroutine fzerox is the main subroutine used for finding the zeros of a com- 
plex function using the Shellman-Morlitt routine as described in Chapter If. Fzerox 
will convert the coordinates of the search rectangle edges into mesh units and proceed 
with the search for sign changes of Im(f)=90 along each edge. ‘The maximum number of 
crossings of Im(=0 allowed with any edge of the search rectangle 1s 100, 

Fzerox allows a maximum number, maxnsy, of mesh squares to lie on any one 
phase curve, Im(M=0. Tf this number is exceeded. the program will reduce the mesh size 


by one-half and start over. If the problem persists after the mesh size has been reduced 


by maxnt times, an error message will be written to the output file and the program 


stops. In the program, maxng is set at four times the number of mesh squares along the 
longest side of the search rectangle. 

Vzerox also checks for modes found on the same phase line. If so, the program 
wil reduce the mesh size and start over. If the problem remains after mesh size has been 
reduced by 2" times. an error message Will be written to output file and the program 
will stop. 

2. Calling statement 
The subroutine ts called by the statement 
Call Fzerox (tleft. tright, tbot, ttop, tmsho, tol, mprint, zeros). 
3. Inputs 
The input variables are 
© tieft : Value of the real part of z at the left edge of search rectangle. 
e «ight » Value of the real part of z at the right edge of search rectangle. 


e thot : Value of the imaginary part of z at the bottom edge of the search 
rectangle. 





e ttop : Value of the imaginary part of 2 at the top edge of the search 


rectangle, 
¢ tmsho : Initial sive of the mesh square. 
® tol : Tolerance to which zeros are to be found. Zeros located closer than 


tol cannot be distinguished. 
¢ mprint > A flag for debugging output. If mprint=0, no debugging will be 
given. Mprint= 1 will acuvate debugging printout. 


4. Outputs 
The output variables are 


® zeros : Array containing the complex zeros of f(z) in the specified 
search rectangle in the complex z-plane. 


e onrz : The actual number of complex zeros of f(z} found. 


§. Calling program element 


® fndmod 
6. Subroutines called 
® findfx 
@ quad 
* nomshs 
7. Common block areas 
®* newmsh 


* tunecom 


N. SUBROUTINE HTGAIN 
1. Description 
This subroutine evaluates the normalized height gain function given by 
Equation 77. 
2. Calling Statement 
The subroutine is called by the statement 
Call htgain (htga. htge. zero, z, m). 
3. Inputs 


The input variables are 


® zero : Modal eigenvalue in g,,-space for mode m. 
e Zz : Height (in meters) above ground at which height gain is to be 
evaluated. 
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m : Mode index. 


4. Outputs 
The output vanables are 
htga : Complex height gain factor in extended arithmetic. 


htge : Real height gain exponent in extended arithmetic. 


5. Calling program element 


modsum 


6. Subroutines called 
xedui 


xcudd 


7. Common block areas 
com? 
pap} 
pap2 
O. SUBROUTINE MODSUM 
1. Description 
This subroutine calculates the field strength relative to free space given by 
Equations 100 and 104. 
2. Calling statement 
The subroutine is called by the statement 
Call modsum (ecms, ems, xm, rngfac, zr. zt, qeigen. ncount). 
3. Inputs 
The input variables are 
xm : Range in meters. 
rngfac : Range factor calculated in the main program. 
zr : Receiver height in meters. 


zt : Transmitter height in meters. 


geigen : Array containing eigenvalues in g,,-space . 


ncount : Counter set in the main program to avoid unnecessary calulations of 
acoefa, acoefe. bcvela and beoefe. If ncount = 1, subroutine 
modsum will call subroutine abcoef. otherwise subroutine 
modsum wil] compute modesum only. 





4. Outputs 


The output variables are 


@ ecms : Coherent field strength relative to free space. 
e cims : Incoherent field strength relative to free space. 


5, Calling program element 


* main 


6. Subroutines called 


abcoef 
* htgain 
® xcadd 


ee addx 


7. Common block areas 


e coml 
¢ com2 
® pap2 


P. SUBROUTINE NOMSHX 
1. Description 
Subroutine nomshx takes the approximate locations of the complex zeros of 
f(z) provided by fzerox and improves the accuracies of the zero locations using New- 
ton-Raphson iteration. 
2. Calling statement 
The subroutine is called by the statement , 
Call nomshx (tol, zeros, nrz). 
3. ‘Inputs 
The input variables are 


e tol : The tolerance to which zeros are to be located. Zeros located closer 
than tol cannot be distinguished. 


¢® zeros : Array containing the approximate locations of the complex zeros of 
{(7). 
e onr : The number of complex zeros stored in array zeros. 
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4. Output 
The output variable is 
e@ zeros > Aus array containing the iterated Values of the complex zeros of f(z). 
5. Calling program element 


e = fzerox 


6. Subroutine called 


e = fdfdtx 
7. Common block area 
© Newmsh 


Q. SUBROUTINE NORME 
1. Description 


This subroutine “normatzes” the complex extended numbers, z = za exp(ze) such that 
exp(—L.0) < max] Re(za)}. | datzad} ] < exp tl.) (342) 


and ze has integer values. 
2. Calling statement 
The program 1s called by the statement 
Call norme (za, ze). 
3. Inputs 
The input variables are 
¢ 7a : Complex amplitude of the extended complex number, z. 


e ze : Real exponent of the extended complex number, Z 


4. Outputs 


The output variables are 


e 7a : Complex amplitude of the “normalized” extended complex number, 
z. 
® ze : Real exponent of the “normalized” extended complex number, z. 


5. Calling program elements 


e xcadd 
e xcedait 
© xcdaig 
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6. Subroutine called 


6 none 


7. Common block area 


® none 


R. SUBROUTINE NORMRE 
1. Description 


This subroutine “normalizes” the real extended number. 7 = za exp(ze), such 


exp(—-1.0) < fray < exp( 41.0) 


and 7e has integer values. 
2. Calling statement 
The subroutine is cailed by the statement 
Call normre (za, ze). 
Inputs 
The input variables are 


: Real amplitude of the real extended number, z. 


: Real exponent of the real extended number, z. 


Outputs 
The output variables are 
: Real amplitude of the “normalized” real extended number, 2. 


: Real exponent of the “normalized” real extended number, z. 


5. Calling program element 


addx 


6. Subroutine called 


none 


7. Common block area 


none 








S. SUBROUTINE QUAD 
1. Description 
This subroutine finds the roots of the quadratic equation of the form given bv 
Equation 2S3. If} ¢ | given in Equation 285 is less than 0.3, the roots are computed with 
Equations 295 and 296. Otherwise Equations 286 and 287 are used. In addition. if 


je! 
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i XS 10-* | subrouune quad returns only I solution given by Equation 293, 
“2, Calling statement 
The subroutine is called by the statement 
Call quad (a, b,c. sol. nrsol. mprint). 
3. Inputs 


The input variables are 


e 4 : The coeflicient of v7 in the quadratic equation. 
e ob : The coefficient of 2x in the quadratic equation. 
ec : The constant term in the quadratic equation. 

q | 


4. Outputs 


The output variables are 


esol > An array containing the real roots of the quadratic equation. 
e nrsol : The number of real roots found. 
* nprint > A flag for debugging. If mprint is not equal to zero, values of the 


coefficients a, b, and ¢ will be written to file.out. 


§. Calling program element 


e =6fzerox 


6. Subroutine called 


e none 


7. Common block area 


e none 


T. SUBROUTINE SEAH20 
I. Description 
This subroutine evaluates the real effective relative dielectric constant and the 
real effective conductivity of scawater as a function of temperature, salinity and fre- 


quency. The program is based on the equations presented in section G of Chapter II. 
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Calling statement 
The subrouune is called bs the statement 
Call seah2o (sigell. epsell. t. s. freq). 
Inputs 
The input variables are 
: Temperature of seawater in degree Celcius. 
: Salinity of seawater in grams of salt per kilogram of seawater. 
: Frequency in Hertz. 
Outputs 
The output variables are 
: Real effective conductivity in siemens per meter. 
> Real effective relative dielectric constant. 
Calling program element 
main 


6. Subroutine called 


none 


7. Common block area 


none 


U. SUBROUTINE SURF 


Description 


ae d;, ; a8 
This subroutine evaluates +, and ——— when surface roughness is included in the 
ad 


model. 

2. Calling statement 

The subroutine is called by the statement 
Call surf (qil, gamma. dgamdq). 

3. Input 

The input variable is 
e gill : Location in the complex g,,-plane where the modal function and its 
derivative are to be evaluated. 

4. Outputs 

The output variables are 


®* gamma : This is ;, given by Equation 180 and 192 for horizontal and vertical 





polarization respectively. 





© dgamdyg + Tins is =— given by Equations 187 and 193 for honzontal and 
cid 
vertical polarizuuon respectively. 
5, Calling program elements 
¢  fetvly 
e  fdfdtx 


6. Subroutines called 


e = function ctunh 


e  funcuen cseche 
7. Common block areas 
¢ com! 
* coms 
* come 


V. SUBROUTINE XNAINEG 
1. Description 
This subroutine evaluates the Airv function, Ai(z7), and its derivative, .31 (2), for 


complex Z in extended complex arithmetic for 





Ss arg{z) «< 0 . (344) 


The subroutine will test if 2* is within the ellipse with foci (—ISS, 3.95) and 
(0.28. —2.11). major axis of 9.0 and minor axis of 64. Hz" is within the ellipse, xaineg 
will invoke xcdait to compute the Airy function and its derivative with a Taylor series. 
Otherwise xaineg invokes xcdaig which uses the Gaussian quadrature method. 
2. Calling statement 
The subroutine is called by the statement 
Call xaineg (z, ata, aie, dala, daic). 
3. Input 
The input variable ts 


e 7 : A value in the complex 7-plane at which the Airy function and its 
derivative are to be computed. 
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Outputs 
aha : Complex amphtude of Az) in extended complex arithmetic. 
aie > Real exponent of Aw) in extended complex anthmetic. 
data > Complex amplitude of Ai‘ (7) in extended complex arithmetic. 


diie : Real exponent of Ai! (4) in extended complex anthmetic. 


8. Calling program element 
xedat 
6. Subroutines called 
xeduit 
xediurg 
Common block area 


none 


SUBROUTINE NAIPOS 
1. Description 
This subroutine evaluates the Airy function. Aw) and its derivauve, Ar(7) in 


extended comples anthmete for 


Q < arg(s) <= : (345) 
£4y 


aq 


The subroutine will test if 7 1s within the ellipse with foci (—1.88. 3.98) and 
(0.28, —2.11), major axis of 9.0 and minor anis of 6.4. 1f 2 is within ellipse, xaipos will 
invoke Acdait to compute the Airy function and its derivative with a Tavlor series. Oth- 
erwise Xaipos invokes xcdaig which uses the Gaussian quadrature method. 

2. Calling statement 
The subroutine ts called by the statement 
Call xaipos (7, ata. aie, data, daic). 
Input 
The input variable is 


> A value in the complex z-plane at which the Airy function and its 
derivative are to be computed. 








4. Outputs 


The output variables are 


ea > Compiles amphtude of Ais) in extended complex arithmetic. 
© ale > Real exponent of Aw) in extended comples aridhiietic. 
@  daia > Complex amphtude of Ai(7) in extended complex arithmetic. 
© date : Real esponent of A(z) in extended complex arithmetic. 


§. Calling program element 


® xvcda 


6. Subroutines cafled 
e xedalt 
® xcdiig 


~” 


7. Common block area 


° none 


X. SUBROUTINE XCDAI 
1. Description 
The subroutine is the driving program in the evaluation of the Aim function, 
Az). and its derivative, AV(7). IF 


GQ < arg(z) < = : (346) 


Acdat invekes AXuipos to compute AZ) and Ai(7) and invokes Xaineg to compute 
eh(ze 74) and (se 773). Pt then uses the connection formula of Airy funcuon in 
Equation 302 and its derivative to compute faze 73) and vfi'{re 24). 

If 


== £ aret) <= 6 . (347) 


Acdai invokes xaipos to compute Adce 274) and .f/'(ze 277). It also invokes xaineg to 
compute Ai7) and Ai(7). It then uses the connection formula to compute -filze ~7"4) 
and -fi(ze 2), 

If 





10s 








or 


—z < arg{z) < --> , (349) 


xedai calls xaipos to compute -fi(ze 7773) and (ce ~774). Tt also calls xaineg to com- 
pute lace 774) and Ar(ze 7°). Tt then uses the connection formula to compute Ai?) 


and A(z). 


if 





(350) 


xcdai calls Xaipos to compute A(z) and Ar(z) and calls Xaineg to compute -filze 77+) 
and (ze? 7*) 2 It then uses the connection formula to conpute faze 777) and 


Ar(ze 74). 


—— < arg(z) < - (351) : 


tof 


xcdai calls xaipos to compute Ai(ce 773) and Ai(ce 75). Tt also calls xaineg to compute 
Anz) and A(z). It then uses the connection formula to compute -fize 774) and 
Aire o 3). 
2. Calling statement 
The subroutine is called by the statement 
Call xcdai (z, aia, aie, daia, daic. aipa, aipe, daipa. daipe. aima., 
aime, daima, daime). 
3. Input 
The input variable is 
e 2 : A value in the complex z-plane at which the Airy function anl its 
derivative are to be computed. 
4. Outputs 
The output variables are 


* ala : Complex amplitude of Ai(z) in extended anthmetic. 
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® aic 

e daia 

e § daic 

¢  aipa 

® aipe 

e daipa 
e daipe 
® aima 
e aime 
e daima 
e daime 





: Real exponent of Ai(z) in extended complex arithmetic. 


: Complex amplitude of Ai‘(z) in extended complex arithmetic. 

: Real exponent of Ai‘(z) in extended complex arithmetic. 

: Complex amplitude of Ai(ze #4) in extended complex arithmetic. 

: Real exponent of Ai(ze-?") in extended complex arithmetic. 

: Complex amplitude of Ai’(ze #"3) in extended complex arithmetic. 
: Real exponent of Ai'(ze #*) in extended complex arithmetic. 

: Complex amplitude of Ai(ze ~?"°) in extended complex arithmetic. 
: Real exponent of Ai(ze ~?"3) in extended complex arithmetic. 

: Complex amplitude of Ai’(ze -?"3) in extended complex arithmetic. 


: Real exponent of Ali’(ze -2"3) in extended complex arithmetic. 
p P 


5. Calling program element 


© = fctvix 


6. Subroutines called. 


® xaipos 


e xaineg 


xcadd 


7. Common block area 


® none 


Y. SUBROUTINE XCDAIG 


1. Description 


The subroutine computes the Airy function, Ai(z), and its derivative, Ai‘(z), in 


extended arithmetic using the Gaussian quadrature method. 


2. Calling statement 


The subroutine is called by the statement 


3. Inputs 


Call xcdaig (2, aia, aic, daia, daic). 


The input variables is 


: A value in the complex z-plane at which the Airy function and its 
derivative are to be computed. 
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4. Outputs 
The output variables are 
ala : Complex amplitude of Ai(z) in extended arithmetic. 
aie : Real exponent of Ai(z) in extended complex arithmetic. 
diua : Complex amplitude of A1(z) in extended complex arithmetic. 


duie : Real exponent of Ai(z) in extended complex arithmetic. 


5. Calling program elements 
AIPOS 


xaineg 


6. Subroutine called 


norme 


7. Common block area 


norme 


SUBROUTINE XCDAIT 
1. Description 

The subroutine computes the Airv function, Ai(7), and its derivative, A1(z), in 

extended arithmetic using a Taylor series expansion, 

2. Calling statement 

The subroutine is called by the statement 

Call xcdait (Z. aia. aie, dala, daie). 

3. Input 

The input variable 1s 

: A value in the complex z-plane at which the Airy function and its 
derivative are to be computed. 

4. Outputs 

The output variables are 


aia : Complex amplitude of Ai(z) in extended arithmetic. 


aic : Real exponent of Ai(z) in extended complex arithmetic. 


gaia : Complex amplitude of Ai‘(z) in extended complex arithmetic. 


daie : Real exponent of Ai‘(z) in extended complex arithmetic. 





5. Calling program elements 
® Xaipos 


® xXaineg 


6. Subroutine called 


® norme 


7. Common block area 


® none 


AA. FUNCTION PROGRAMS 
1. Function casin 
Function casin evaluates the complex arcsine of a complex number, z. with the 
equation 
sit fe) = —fdn fist)? Se ey) 5 (382) 
The function is called by the statement 
casin (2). 
The function is used in the main program to convert the eigenvalues from 
G,,-plane to @-plane . 
2. Function csech2 
Function csech2 evaluates the square of the hyperbolic secant of a complex 


number, z. with the equation 
4 en. 
sech* (=) = oF (353) 


The function is called by the statement 


csech2 (2). 








: . : : Ca Ch Ca Ch 
It is used in subroutine surf for evaluating — ae gage ANG Se 
CHa CF, CGin CF 


(see Section F, Chapter I). 
3. Function ctanh 
Function ctanh evaluates the hyperbolic tangent of a complex number, z, with 


the equation 
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tanh (z) = ——“ 
| oes 
: : atm ee : 
for | Re(z) | > 0 as (354) 
or | Jaz) | > | 
lo .2 
or 
1 3 2 6s [?. 3 6: 9 
tanh (zs) =z —- me Ts 7 7 3p? + Was Fo 
| 
Az < = Pee 
for | Re(z) | < ra (345) 
] 
or | dmc) | < — 
[<2 


The function is called by the statement 
ctanh (2). 
It is used in subroutine surf for evaluating a. b, @ b and their derivatives. 
4. Function ibstrip 

Function ibstrip finds the length of a character string by removing all trailing 
spaces from the end of the string and returns the number of characters in the stripped 
string. 

The function is called by the statement 

ibstrip (inst, outst). 

Where inst is the input string and outst is the stripped string. 


The function is used in the main program to check the length of the filename. 


AB. MLAYER SUPPORTING PROGRAMS 
MLAYER has several supporting programs that helps in the running of the pro- 
gram. They are 
e mlaver.mak 
e wvgclcan.mak 
© wvgstrip.mak 


e wvegclean.c 
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e  wvegstrip.for 


e = wvestdin.for 


1. Miayer.mak 
This is the makefile for the creation of mla20 and mlaS80 programs. Mla20 uses 
subroutine fndmod20.for to set up the rectangular search region Whereas mlaSO uses 
fndmods0.for. 
Miaver.mak compiles and links all program elements and creates mla20.exe and 
miaS0.exe for the execution of MLAYER program with either mla20 or mlaS0o. 
To compile and link MLAYER, type 
nmiuke miaver.mak 
Nmake will Compare the modification dates of the target files (e.g.. object and executable 
files) with those of the dependent files (e.g., source files), If any of the dependent files 
has been changed recently, 1e., the modification dates of the dependent files are more 
recent than the target files. nmake will execute the commands in miaver.mak to update 
the target files. Nmake updates only those outdated target files. 
To execute MLAYER program. type 
mla20 < infile 
to execute the mla20 version. Infile is the name of the input file which can be either filein 
or filein.eig. If filein is used, MLAYER will search for all modes of propagation. If fi- 
lein.eig is used. MLAYER will compute the modsum and path loss without a search for 
the modes. 
2. Wrvygclean.mak 
This 1s the makefile for the creation of wvgclean.exe 
3. Wvygstrip.mak 
This is the makefile for the creation of wvgstrip.exe. 
4. Wvgclean.c 
Wvgclean.c adds comment or descriptors to each datum line in filein and fi- 
lein.eig. This will help in editing. 
To execute the program, type 
wvgclean < filein.cig > temp 
ren temp filein.eig 
5. Wygstrip.for 


Wvestrip.for removes comments or descriptors that are added by wvgclean.c. 


To execute the program, type 











wvestrip < fileineig > temp 
ren temp filein.cig 
6. Wvegstdin.for 
Wvestdin.for reads input data from filein or filein-eig on execution of the 
MLAYER program. 


VI. DISCUSSIONS AND RECOMMENDATIONS 


A. ASSESSMENT OF MLAYER 

MLAYER was developed with the intention of providing a program which is ca- 
pable of locating all propagating modes with attenuation rates below a specified value. 
Jlowever, the price of this capubilitv is the extremely Jong execution time. For instance, 

a 2-meter evaporation duct takes approximately 3 hours and 20 minutes (on an IBM 
~ PS 2 model SU with an Intel $0386 processor at 16 MIIz) to locate 9 modes with atten- 
uation rates less than S decibels per kilometer. A }4-meter evaporation duct takes about 
6 hours to locate 94 modes with attenuation rates below 2.1 decibels per kilometer. 
Samples of output data are found in Appendix B. 

Calculations of field strength and path Joss with MLAYLR were found to agree 
favourably with experimental measurements [Ref 17 and IS}. Typically, the model un- 
Jerestimates the measured valucs by about 10 dB. The sources of discrepancies were 
probably due to the validity of the surface roughness model and the assumption of lat- 
eral homogeneity of the medium. 

In the following sections, proposed areas for euhancement of the capabilities of 
MLAYER are discussed. 


B. SURFACE ROUGHNESS MODEL 

The surface roughness model used in MLAYER is based on the Kirchhoff approx- 
imation Which ts a single scattering theory. This model is only valid for gently undulating 
surfaces where shadowing of incident and scattered fields are negligible. When shadow- 
ing cannot be neglected, the Kirchhoff approximation tends to overpredict the scattered 
energy Which leads to erroneous results. 

Shadowing occurs when incoming or outgoing ravs are blocked by the rough sur- 
face. This leaves some area of the surface in shadow. This effect is associated with mul- 
tiple scattering. Therefore. to account for shadowing. multiple scattering interactions 
must be included in the model. 

A more accurate model called the phase perturbation technique was introduced by 
Winebrenner and Ishimaru [Ref. 19]. The phase perturbation theory incorporates all the 
multiple scatterings required for shadowing of the incident field. The theory is based on 
the extinction theorem which is a consequence of Green's theorem [Ref. 20]. In Refer- 


ence 19, numerical results of the phase perturbation bistatic scattering cross section were 
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compared with exact numerical results. In the region examined, it was found that the 
results of the phase perturbation technique agree with the exact results except at low 
grazing angles. 

The phase perturbation technique is valid for a wider range of applications. It was 
shown in Reference 19 that the reflection and backscattering coefficients reduce to that 
of the Kirchoff approximation when /// 2 1, ie.. for gently undulating surfaces, where 


/ is the length and « is the wavelength of the incident wave. 


C. LATERAL HOMOGENEITY OF REFRACTIVITY PROFILE 

In MLAYER, the refractive index is assumed to vary only with height. Although 
this assumption appears to be adequate most of the time [Ref. 17]. there are incidents 
Where the validity of lateral homogeneity was cited as a possible reason for discrepancies 
between observations and predicted results [Ref 18 and 21]. 

For refractivitv profile that varies along the path. the laterally inhomogencous 
structures can be approximated by several homogencous sections. The exchange of en- 
ergy between the propagating modes at the boundary of two sections are analyzed by 


mode-conversion techniques [Ref. 22]. 


D. FIELD CONTRIBUTION FROM THE SOURCE 
In the formulation of the propagation model, MLAYER neglects the contribution 


from the direct wave of the source. This approximation is good for ranges near or be- 


vond the horizon where the major field contributions are from the reflected waves. For 
ranges Well within the horizon, the contribution from the direct wave is significant and 
should not be ignored. Hence. for wave propagation in the interference region, the co- 
efficients, 4, and B, of Equations 31 and 32 must be obtained from the solutions of 
Equation 46. At the same time, the surface wave field contribution from the integral of 


the contour, C, , in Equation 73 may become significant. 








VI. CONCLUSION 


MLAYER is a useful research tool for conducting case studies over a large dynamic 
range of frequencies. Calculation of field strength and path loss are in agreement with 
observations for frequencies from 63 Mhz to 94 GHz [Ref. 17 and 1S]. It can serve as a 
yardstick against which results of other quicker but less accurate models such as [REPS 
[Ref. 23] and LREPS [Ref. 24] are compared. In addition. the multiple (as opposed to 
trilinear) piecewise linear refractivity profile makes it possible to model the simultaneous 
occurrence of an evaporation duct and elevated ducts. 

Finally, modifications of MLAYER to incorporate conditions of horizontal heter- 
ogencity of refractivity, effects of source contribution and an improved mede!ling of 


surface roughness as discussed in Chapter VI are recommended to fully exploit its ca- 


pabilities. These are oppértunities opened for further studies. 











APPENDIX A. REPRESENTATION OF &, (¢) AND k,(q) IN TERMS OF 
MODIFIED HANKEL FUNCTION OF ORDER ONE-THIRD 


The Airv function is related to the modified Bessel function of the third kind of or- 


der one-third by the following equation [Ref 4]: 


Ai = + (Ey Kia (Fe (356) 
where 
Ky3(g) = fe" ha (ge) (357) 
or 
Ky lq) = i ey Oger) (358) 


In Equations 357 and 358, //,,;%(.) and //,,®(.) are Hankel functions of the first and 


acn7 


second kind of order one-third respectively. Substituting Equation 357 into Equation 336 


leads to 
¢ aa -” jm! . ' poets, 
Ai(g) = = (4 ty eeh in (are) 
(339) 
} YJ \12 z'6 di) 2 1-999 ooo 
am aL ead 
Thus, 
Ai( —ge?*?) fs (4)! (e -jr yi (e273 IP 6 ini6 
()f 2 3/2 35/2 je jaf? 
x My, 79 € ere (360) 
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From Equanon 142, 
Kiss =p TA lt See: (361) 


Subsutuung Equation 360 into Equation 301 gives 


io f \1'2 a) oa 
Ay (q) = y2yh (= ye hs" (> _) 
sgh? 8 ar 362 
= (+ pee Hf as (+ q *) (362) 
= hh, (qq) 


Where A, {y) is the modified Hankel function of the first Kind of order one-Uurd and is 


defined as [Ref. 4] 


> mo \ae g : 2 
Ay(g) = ‘ey gaye Tn 7’) 103) 


To derive an expression for A. (gy), the Airy function, -f(— gh can be expressed as 
[Ref. 4] 


ss oe Lf @\ve2| iste il an ee Gee oy poe : 
li -Qg= a cae) E ; Hy: At q ' ee nig (eee )] (3604) 
From Equation 143, 


kstg) = 2(12)"* ef" Ail —q) . (368) 


Substituting Equation 364 into 36S leads to 


117 





9 an { jo! 1 % "y Do) ry 

= eee B72 \1/3 | d5/3 y (i) = gory aT 4 2 422 
y f 1/3 oy { cn 3 (360) 
i. SOO 


=h;(q) - ese hig), 


Where fy (g) is the modified Hankel function of the second kind of order one-third and 


is defined as [Ref. 3] 
a) Nts sf oD . 
(+ q a cue a r*) (307) 


For large | y|. the Hankel function can be approximated by the leading term of its as- 


hy (7) 


Mf 


ymptouc expansion. That is |Ref. 4]. 


Hae) ~ aye ap[i(y- 2-4) | 





(368) 
—m < arg(y) < 2n 
and 
my ~ (Sr) eo[-iv-2-4)] 
173 mq P ae 6 4. (369) 


LIS 














; ees ; > . S25 a 
ai-6 (-i2 Als ? we ee. a . (370) 
raeaal baller g exp A Se )] ; 


< arg{q) < — + 





To obtain the asymptotic expansion for kyq) . it is more convenient to express A-fg) In 


terms of J/,,%(.). From Equations 356 and 388, 


1374) 
Tf GN ays mf 2 32 a3 
=> ( ~ ) Ce uw ( 4 wo 
co by : = 
Thus, 
Df G Nin jad mjete yf 2 22 [niet oni 
A —g) = -+(+) Tee = oF (Sg ( ( ) 
(aon) 
! ( q \ 2 jen Git 2 se ae) 
= 5 3 Lies : gy € ‘ 
Substitute Equation 372 into 365 gives 
2 3/2 \172 af 2 3/20 ajas a-4 
ky(q) = > (+ q yy ae (373) 


From Equation 369 and 373, the asymptotic expansion for 4; (q) Is 
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APPENDIX B. 


SAMPLE PRINTOUT OF INPUT AND OUTPUT FILES 


vert Sample printout for a 2 m evaporation duct input filers 


9ghz02m 

0 

9600. 00 
0. 000000 
1. 000000 
0. 000000 
5.000000 
15. 000000 
35. 000000 
. 000000 
. 000000 
. 000000 
. OGQ000 
»25ue850 
48. 000000 
0. 000000 
1. 000000 


ODOooor 


46. 0000000 


7. 0000000 
2.009000000 
18.5960 
9.250000 
3. 000000 
C. 000000 
16 

0. 000000 
0. 000000 
0.000000 
0. 020000 
-0. 400000 
0. 000000 
0. 040000 
-0. 550000 
0. 000000 
0.079090 
-0. 700000 
0.000000 
0. 158000 
-0. 850000 
0. 000000 
0.251000 
-0. 940000 
0. 000000 * 
0.631000 - 
-1. 100000 
0. 000000 
1. 259000 





FQMZIN 
DELFQ 
NFREQ 
MPOL 
ALOSS 
SEATMP 
SEASLT 
IFLAGE 
ALRTMP 
RH 
WGPM3 
RMSBHT 
aTINIT 
DELZT 
NZT 
ZRINIT 
DELZR 
NZR 
XINIT 
DELX 

NX 

ZREF 
NZLAYR 
zi[ 0] 
zin[ 0] 
zigab{ 0] 
zif{ 1} 
zin[ 1] 
zigab[ 1] 
Zi{ 2] 
zim{ 2] 
zigab| 2] 
zi[ 3] 
zim{ 3] 
zigab{ 3] 
zi[ 4] 
Zim| 4] 
zigab[ 4] 
zi{ 5} 
zim[ 5) 
zigab[ 5] 
zil[ 6] 
zim{ 6} 
zigab[ 6] 
zi[ 7] 


-1. 170000 
0. 000000 
2. 000000 
-1. 180000 
0. 000000 
3.981000 
-1. 090000 
Q. 000000 
5.012000 
-1. 010000 
0. 000000 
7.943000 
-0. 750000 
0. 000000 
12. 589000 
-0. 270000 
0. 000000 
19. 953000 
0. 550000 
0. 000000 
25. 119000 
1. 140000 
0. 000000 
31. 623000 
1. 900060 
0. 000000 
50. 000000 
4.069000 
0. 000000 


dededededededed S amp le pr intout for 9 gh z02.e i givievediederere 


9ghz02m 
1 


9600. 000 

0. 000000 

1 

0 

5. 000000 

15. 000000 
35. 000000 
1 

0. 000000 

0. 000000 

0. 000000 

0. 250000 

48. 000000 
0. 000000 


1 

46. 060000 
7.000000 
2 

18. 500000 
9. 250000 
3 


zinf{ 7] 
zigab[ 7] 
zi| 8] 
zim{ 8] 
zigab[ 8] 
2i[ 9] 
zim[ 9] 
zigabl 9] 
zi[ 10] 
zim[ 10] 
Zigab[ 10} 
zi[ 11] 
zim[ 11} 
zigab[ 11} 
zi[ 12] 
zim{ 12} 
zigab[ 12] 
zi[ 13} 
zim{ 13] 
zigab[ 13] 
zif 14] 
Zim[ 14} 
zigab[ 14] 
zi[ 15] 
zim[ 15] 
zigab[ 15] 
zi[ 16] 
zim[ 16] 
zigab[ 16] 


FQMZIN 
DELFQ 
NFREQ 
MPOL 
ALOSS 
SEATMP 
SEASLT 
IFLAGB 
AIRTMP 
RH 
WGPN3 
RMSBHT 
ZTINIT 
DELZT 
N2T 
ZRINIT 
DELZR 
NZR 
XINIT 
DELX 
NX 





0. 000000 
16 
0. 000000 
0. 000000 
0. 000000 
0. 020000 
-0. 400000 
0.000000 
0. 040000 
-0. 550000 
0. 000000 
0. 079000 
-0. 700000 
0. 000000 
0. 158000 
-0. 850000 
6. 000000 
- 251060 
-0. 940000 
0. 000000 
0. 631000 
-1. 100000 
0. 000000 
1. 259000 
-1. 170000 
0. 000000 
2.000000 
-1. 180000 
0. 000000 
3. 981000 
-1. 096900 
0.000000 
5.012060 
-1.010000 
0.000000 
7.943000 
-0. 750000 
0. 000000 
12. 589000 
-0. 270000 
0. 000000 
19, 953000 
0. 550060 
0. 000000 
25. 119000 
1. 140000 
0. 000000 
31. 623000 
1. 900000 
0. 000000 
50. 000000 
4.069000 
0. 000000 
9 


ZREF 
NZLAYR 
zi[ 0] 
zim|[ 0] 
zigab[ 0] 
zi[1] 


zin{ 5] 
Zigab[ 5] 
zi[ é] 
zinm[ 6] 
Zigabf 6] 
Z2i[ 7] 
zim[ 7] 
zZigab[ 7] 
zil[ 8} 
Zim{ 8} 
Zigab[ 8] 
zi{ 9} 
Zim{[ 9] 
zigab[ 9] 
2Zi[ 10} 
zZim[ 10] 
zigab[ 10] 
zi{ 11] 
zim[ 11] 
Zigab[( 11] 
zi[ 12] 
zimf{ 12] 
zigab[ 12] 
zi[ 13] 
Zim[ 13] 
zigab[ 13] 
zi[ 14] 
zim{ 14] 
zigab[ 14] 
zi[ 15] 
Zim[ 15] 
zigab[ 15] 
zil[ 16} 
zim[ 16] 
zigab[ 16] 
nrmode 


(-1. 269556983588969E-001,1.637613838717031E-001) 
(-9. 291194300692691E-002,1. 251243717867209E-001) 
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oe 
Low od 





ee et 








(-4. 143216670274859E-002,8. 947793427709196E-002) 


(8. 240706115401963E-003,6. 
(5. 551160395245244E-002,5. 


a ae a a 
Ot oe 


1 


frequency = 


- 6826113485045449E-002,6. 
. 394420773218969E-001,9. 
. 102775842495503E-001,1. 
. §864589183106626E-001, 1. 


9600. 0000 mhz 


horizontal polarization 


maximum mode attenuation rate = 


seawater temperature = 


seawater salinity = 
dielectric constant of seawater 
conductivity of seawater = 


rms surface bump height 


iflagb = 


air temperature 
relative himidity = 


917378956118456E-002) 
324818113693566E-002) 
708823745503940E-002) 
725689299213881E-002) 
357828223763695E-001) 
658090146540282E-001) 


Se ere 
WmonNInAM Ew 
ee ss a 


5.0000 db/km 


.0000 degrees celsius 
.0000 percent 


liquid water concentration in air = 


(meters) 


. 0000 
. 0200 
. 0460 
0790 
. 1580 
2510 
. 6310 
1. 2590 
2.0000 
3.9810 
5.0120 
7.9430 
12. 5890 
19. 9530 
25.1190 
31. 6230 
50. 0000 


m(z) 
(m-units) 


- 0000 
- 4000 
- 5500 
. 7000 
- 8500 
- 9400 
. 1000 
. 1700 
- 1800 
. 0900 
. 0100 
. 7500 
. 2700 
- 5500 
. 1400 
- 9000 
- 0690 


. 0000 m-units 
.0000 db/km 


. 00000D+00 
. 00000D+00 
- OOOODOD+00 
. 00000D+00 
. OVOO0OD+00 
. 00000D+00 
. OOOOOD+00 
- 00000D+00 
. 00000D+00 
- OOODOD+00 
. OOOO0OD+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. OOOOOD+00 
. OOODOD+00 
. O0000L+00 
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-20. 
-7. 
-3. 
1. 

9677 

~4211 

~12115 

. 0135 

0454 

0776 

. 0887 

. 1033 

1114 

~ 1142 

- 1169 

. 1180 


15.0000 degrees celsius 
35.0000 grams salt/kg seawater 
= 80. 8869 

4.6400 si/m 
. 2500 meters 


tropospheric modified refractivity profile 


gab(z) 
(db/km) 


dm/dz 
(m-units/meter) 


0000 
5000 
8462 
8987 








.0000 grams/meter***3 


d( gab) /dz 
((db/km) /meter) 


. 90000D+00 
. 0OOVOD+00 
. 00000D+00 
- 90000D+00 
- 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
- 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. QOOOOD+00 





te ateatan’, 


tmesh= 1.67984D-03 


zeros found in expanded 


tleft= 0. 0OD+00 

ttop = 1. 68D-01 
geigen( 1) = 
geigen( 2) = 
qeigen( 3) = 
qeigen( 4) = 
qeigen( 5) = 


zeros found in expanded 


tleft= -2.69D-01 

ttop = 1. 68D-01 
geigen( 1) = 
qgeigen( 2) = 
qeigen( 3) = 


zeros found in expanded 
tleft= -5. 38D-01 
ttop = 1. 68D-01 


zeros found in expanded 
tleft= -8. 06D-01 
ttop = 1. 68D-01 


zeros found in expanded 


tleft= 2. 69D-01 
ttop = 1. 68D-01 
qeigen( 1) = 


zeros found in expanded 
tleft= 5. 38D-01 
ttop = 1. 68D-01 





Start search for modal eigenvalues wereaerede 


tol= 1. OOO00D-04 


search rectangle defined by 


tright= 2.69D-01 

tbot = -1.00D-04 
8.24071D-03 6.91738D-02 
5.55116D-02 5.32482D-02 
2.10278D-01 1.35783D-01 
1.39442D-01 9.72569D-02 
7.68261D-02 6. 70882D-02 


search rectangle defined by 
tright= 0. OOD+00 
tbot = -1. 00D-04 


-1. 26956D-01 
-9,29119D-02 
-4.14322D-02 


1.63761D-01 
1.25124N-01 
8. 94779D-02 


search rectangle defined hy 
tright= -2.69D-01 
tbot = ~1. 00D-04 


search rectangle defined by 
tright= -5. 38D-01 
tbot = -1. 00D-04 


search rectangle defined by 
tright= 5. 38D-01 
tbot = -1. 0OD-04 


2. 86459D-01 1. 65809D-01 


search rect*ngle defined by 
tright= 8. 06D-01 
tbot = -1. OOD-04 


~~ 
ts 
a 





zeros tounu in expanded search rectangle defined by 


tleft= 8. 06D-01 tright= 1. O8D+00 
ttop = 1. 68D-01 thot = -1. 00D-04 
yerererere modal eigenvalues nrmode= = 9 ¥iritic 
mode eigenvalue theta atnu db/km 
1 -1.26956D-01 1.63761D-01 1.16912D-03 2.38566D-03 4.8743D+00 
2 -9.29119D-02 1.25124D-01 1.03533D-03 2.05835D-03 3.7243D+00 
3 -4.14322D-02 8.94779D-02 9.86784N-04 1.54436D-03 2. 6633D+00 
4 8.24071D-03 6.91738D-02 1.15188D-03 1.02280D-03 2. 0589D+00 
5 5.55116D-02 5.32482D-02 1.50185D-03 6. 03859D-04 1. 5849D+00 
6 7.68261D-02 6.70882D-02 1.74517D-03 6.54734D-04 1. 9969D+09 
7 1.39442D-01 9.72569D-02 2.29575D-03 7. 21529D-04 2. 8948D+00 
8 2.10278D-01 1.35783D-01 2.80081D-03 8.25696D-04 4.0415)+00 
9 2.86459D-01 1.65809D-01 3.24286D-03 8&.70843D-04 4.9353D+00 
sevedededeterede S amp le pr intout for 9 gh 202m. P Le redsrriedevededs 
3 1 2 
frequency = 9600.0000 mhz 
18.5 48.0 46.0 52.90 61.21 84.54 76. 23 
18.5 48.0 53.0 69.58 73. 08 67. 86 64. 36 
27.8 48.0 46.0 2.17 38.48 138. 80 102. 48 
27.8 48.0 53.0 19.57 43.39 121. 39 97.57 
37.0 48.0 46.0 3.41 22.15 140. 05 121.31 
37.0 48.0 53.0 2.69 25.91 140.77 117.55 


ieee Sample printout for a 14 m evaporation duct input filerrereredrerer 


9¢ghz14m 

0 

9600. 00 FQMZIN 
0.000000 DELFQ 
1 NFREQ 
0 MPOL 
2.100000 ALOSS 
15.000000 SEATNP 
35.000000 SEASLT 
1 IFLAGB 
0.000000 AIRTMP 
0.000000 RH 
0.000000 WGPM3 
0.250000  RNSBHT 
25.000000 ZTINiT 
0.000000 DELZT 
1 NZT 
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DUanwnou 





3. 0000009 
7.0000000 


0. 000000 
0. 040000 
-3. 900000 
0. 000000 
0. 100000 
-5. 340000 
0. 000000 
Q. 200000 
-6. 400000 
0. 000000 
0. 398000 
-7. 460000 
0. 0000N0 
0.794090 
-8. 490000 
0. 000000 
1. 585000 
-9,.470000 
0. 000000 
3. 162000 


-10. 350000 


0. 000000 
6. 310000 


-11. 040000 


0. 000000 
12. 589000 


-11. 320000 


0. 000000 
14. 000000 


-11. 320000 


0. 000000 
25. 119000 


-10. 870000 


0. 000000 
39. 811000 
-9. 750000 
0. 000000 
50. 119000 
-8. 820000 
0. 000000 
63. 096000 
-7. 560000 
0.000000 
79. 433000 
-5. 880000 
0. 000000 


ZRINiT 
DELZR 
NZR 
AINIT 
DELX 
NX 


zim[ 2] 
zigab[ 2] 
zi[ 3] 
zim{ 3} 
zigab[ 3) 
zi[ 4] 
zim 4] 
zigab{ 4]: 
zif 5] 
Zim{ 5) 
zigab[ 5] 
z2if{ 6] 
zin|[ 6] 
zigabf 6] 
zif 7} 
zim 7] 
Zigab[ 7] 
zi{ 8] 
zim{ 8] 
zigab| 8] 
2i[ 9] 
zim| 9] 
Zigab[ 9} 
zi[ 10} 
zim[ 10] 
zigab| 10] 
zif{ il] 
zim{ 11] 
zigabf 11] 
zif[ 12) 
zim[ 12] 
Zigab[ 12] 
zi[ 13] 
zim{ 13] 
Zigab[ 13} 
zi[ 14] 
zim[ 14] 
zigab{ 14] 
z2i[ 15] 
zim[ 15) 
zigab[ 15] 
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100. 000000 
-3. 670000 
0. 000000 
125. 655006 
-0. 8000C0 
0. 000000 
158. 489000 
2.920000 
0. 000000 
199. 526000 
7.690000 
0. 000000 
209. 526000 
8. 870000 
0. 000000 


dervedesetevededeS amp le pr intout for 9ghz 14m.e Lg torrerertaeresy 


9ghz14m 

1 

9600. 000 
0. 000000 
1 

0 

2. 100000 
15. 000000 
35. 000000 
1 

0. 000000 
0. 000000 
0. 000000 
0. 250000 
25. 000000 
0. 000000 
1 

3. 000000 
7.000000 
2 

18. 500000 
9. 250000 
3 

0. 000000 
20 

0. 000000 
0. 000000 
0. 000000 
0. 040000 
-3. 900000 
0. 000000 
0. 100000 
-5. 340000 
9. 000000 
0. 200000 
-6. 400000 
0. 000000 





zi[ 16] 
Zin| 16} 
zigab[ 16] 
24, 17) 
Zim[ 17] 
Zigab{ 17] 
zil 18] 
zim[ 18] 
zigab[ 18] 
zi[ 19} 
zim[ 19] 
zigab[ 19] 
zif 20} 
zim[ 20] 
Zigab[ 20] 


FQMZIN 
DELFQ 
NFREQ 
MPOL 
ALOSS 
SEATMP 
SEASLT 
IFLAGB 
AIRTMP 
RH 
WGPM3 
RNSBHT 
ZTINIT 
DELZT 
N2T 
ZRINIT 
DELZR 
N2R 
XINIT 
DELX 

NX 

ZREF 
NZLAYR 
zil 0) 
zin|[ 0] 
zigab[ 0] 
zi[ 1) 
zin{ 1} 
zigab[ 1] 
zi[ 2] 
Zin 2] 
zigab[ 2] 
zi[ 3) 
zim{ 3] 
Zigab[ 3] 
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0. 398000 
-7. 460000 
0. OO000O 
G. 754090 
-8. 490000 
9. 000000 
1.585000 
-9. 470000 
0. 000000 
3. 162000 
~10. 350000 
0. 000090 
6.310000 
-11. 040000 
0. 000000 
12. 589000 
-11. 320000 
0. 000000 
14. 000000 
-11. 3390000 
0. 000600 
25. 119000 
-10. 870000 
0. 000000 
39. 811000 
-9. 750060 
0. 000000 
59. 119000 
-8. 620000 
0. 000000 
63. 096000 
-7. 560000 
0. 000000 
79. 433000 
-5. 880000 
0. 000000 
100. 000000 
-3. 670000 
0. 000000 
125. 893000 
-0. 800000 
0. 000000 
158. 489000 
2.920000 
0. 000000 
199, 526000 
7.690000 
0. 000000 
209. 525000 
8. 870000 


zi[4] 
zim] 4] 
zigab| 4] 


z2i[ 7] 
zim 7] 
zigab[ 7] 
zif 8} 
zim{ 8] 


zi{ 13] 
zim[ 13] 
zigab[ 13] 
zi[ 14} 
zim{ 14] 
zigabf 14] 
zZif 15} 
zimf[ 15] 
zigab[{ 15] 
zi[ 16} 
zim{ 16] 
Zigab| 16] 
zi[ 17] 
Zim[ 17} 
zigab[ 17] 
zi[ 18) 
zim[ 18] 
zigab[ 18] 
zif{ 19] 
Zim[ 19] 
zigah[ 19] 
zi[ 20] 
Zim[ 20] 
zigab|[ 20] 
nrmode 


4 
(-1. 141392174382303E-001,2. 381981532549529E-002) 
(-9. 791475196101274E-002,2. 146753762359050E-002) 
(-8. 353759867176463E-002,1. 997153241435371E-002) 
(-7. 039783167095078E-002, 1. 889606906408722E-002) 
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FO LO LO GO LO LO LO LO LL LO LON LO LD LOL LO LO LON OO LO BOE LO LO LE LOE LOLOL OL LO LO LOLOL LO LOL LO LO LO LO LO LO PO LOR LON LOLA LON LOL LOO 


WWM NH NY DD DYN NM NY NHP NHN NH NHNYD NKR RR RRB Ree HS Pe RRP OW OMNYYDU PF WWHH Wd 1 


Wie poi Ey 


-551321843705835E-003,1. 
- 275397200580210E-002,1. 
- 182707485088627E-002,1. 
. 011306010124748E-002,1. 
- 842803342065884E-002, 1. 
- 694372655547260E-002,1. 
-485441606530146E-002, 1. 
. 227012336327741E-002,1. 
- 029778501629384E-002,1. 
- §603056175576653E-002, 1. 
-555171301174155E-002,1. 


. 226312048960586E-002, 


- 940292680300993E-002,1. 
. 068485002790371E-001,1. 
- 138252208858626E-001,1. 
. 207736300828756E-001,1. 
. 270094738812559E-001,1. 
. 334153343077380E-001 , 1. 
. 403978164901863E-001,1. 
-467611452238766E-001,1. 
- 534013614885507E-001, 1. 
- 596149004805981E-001,1. 
- 652780520579286E-001,1. 
- 709514646617846E-001,1. 
. 772646284812436E-001, 1. 
- 836482895835790E-001, 1. 
- 891728805710848E-001,1. 
-951501729145819E-001,1. 
- 006781595657978E-001,1. 
-059566903740310E-001,1. 
- 106263141043572E-001,9. 
-157518751516274E-001,8. 
. 213117492755677E-001,9. 
. 273294238073336E-001,8. 
. 303987147631242E-001,1. 
- 318960198770044E-001,8. 
- 348881611845628E-001,7. 
-402844396949530E-001,8. 
-455147471247609E-001,9. 
-499562362348160E-001,1. 
- 545987673232816E-001, 1. 
» 600542439477720E-001,1. 
- 662691013974199E-001,1. 
- 726418617644147E-001,1. 
- 787900254603262E-001, 1. 
-849476164544477E-001,1. 
- 913747922726816E-001,1. 
- 980151522337324E-001, 1. 
. 042875049575047E-001,1. 
. 113767357254904E-001,1. 


. §14228837496516E-002,1. 806765450625147F-002) 
- 656628790733621E-002,1. 728118397954074E-002) 
- 585330645452212E-002,1. 66356316120551GE-002) 
$27742000221239'5-002,1. 6283482755016 756-0625 
-546901563155011E-002,1.571655186158130E-002) 
-541956071147925E-003, 1. 540239335815212E-002) 


518851207972415E-002) 
460556299546134E-002) 
455854335720260E-002) 
427613140712134E-002) 
365892306660126E-002) 
365320803509313E-002) 
365090923924493E-002) 
336276100865187E-002) 
282089227646639E-002) 
285572220605075E-002) 
287860792288040E-002) 


- 278978888881959E-002) 


210154187763285E-002) 
208486102139816E-002) 
195444121483925F-002) 
220952877564303E~-002) 
222384497715726E-002) 
139158734459426E-002) 
116571468633296E-002) 
129961336802001E-G02) 
116543666415548E-002) 
136779885138012E-002) 
153156478829780E-002) 
066116492781759E-002) 
023671860551962E-002) 
049763117145163E-002) 
025773869961451E-002) 
007563477610222E-002) 
041321459422345E-002) 
043172189566190E-002) 
926428640075490E-003) 
79§215719693539E-003) 
015024671611194E-003) 
321071087819404E-003) 
185554142465963E-003) 
464799970637166E-003) 
882321847941567E-003) 
719487457459410E-003) 
692237120376304E-003) 
010888147000029E-002) 
010149239614689E-002) 
049907356151861E-002) 
044438363485770E-002) 
092098859083637E-002) 
175196099261803E-002) 
227969492073185E-002) 
246149094591169E-002) 
329730594318235E-002) 
319657577519411E-002) 
314377702678184E-002) 


Or OWON TUR Wh | DUO MAIHDUE WOH OOM BUE Who eH Ore 
as aa a nd rss a wa ea ss al Ol ed ld 


BE WWW WWW WWW WNNMRMNNNNNYN HHP R HRP RPP RHE RO ONIN 


i ae Oe ee 
DD Ur & Ww 


(amp Lauen Leen Lanene Looe) 
oS, ie ae ae 
RP Ow Os! 
eS 


[52] 
{ 53] 


wn 
= 


[55] 
[56] 


uuu 
sO CO st 


[ 60] 








AAATVUANMNAANUAME EEE HERR PR RWWWWWOWWWWH Ww 


- 187819780335167E-001,1. 333680648010814E-002) 
. 264578839391167E-001,1. 384739430048459E-002) 


-41uy4ic000d4528E-001,1. 490152636527356E-002) 
»487677270086091E-001,1. 524482539938137E-002) 
- 563811238953353E-001, 1. 551353021082856F -002) 
- 643646265942290E-001,1.538942161737147E-002) 


- §14931085621074E-001, 1. 618605966683290E-002) 
- §98891425944480E-001,1. 689462812245438E-002) 
- 980319033943261E-001,1. 714047 724320524E-002) 
- 067283853697590E-001,1. 729643131530817E-002) 
- 154949433609382E-001,1. 751352926449171E-002) 
- 246400421526860E-001,1. 750551050950182E-002) 
- 341831113214738E-001,1. 789042401852184E-002) 
-435147531419498E-001,1. 866150403101799E-002) 
- 52528998194615355-001,1. 909787541066863E-002) 
- 618630098407780E-001, 1. 919260738082309E-002) 
-716161770152394F-001,1. 927318684881558E-002) 
- 815020432944205E-001,1. 952391397029549E-002) 
- 917304699196995E-001,1. 972766554431011E-002) 
- 020225407049786E-001,2. 044077417123271E-002) 
- 118252276793764E-001,2. 102681780869346E-002) 
. 218948685680153E-001,2. 102179462275660E-002) 
- 324815199710420E-001, 2. 111602176832850E-002) 
- 432604196740195E-001,2. 136847129075478E-002) 
-541811165125796E-001,2. 172276354461229E-002) 
- 651832495410506E-001,2. 227384315438282E-002) 


- 961376451733623E-001,2. 292879795570507F-002) 
. 097267273995461E-001,2. 326416356305734E-002) 
. 213712082142159F-001, 2. 370559660462316E-002) 
- 330414592077247E-001,2. 426761624643026E-002) 





338112597216490E-001,1.464266747417008E-002) 


729775181063331E-001,1. 559989625580370E-002) 


A ee ed ed Ne Ld Cd Md ee 


758666247777175E-001, 2. 293313938401881E-002) 
666533508029618F-001,2. 286336319718634F-002) 


ne en mn ney en nn em pti pe ene ny eR py ee py FY ts GI GE GS Prey gy PE 


WOUOWUWOW WHOM OHWMOM WM NISSEN SINE SI DRHDADDA DAA 
Pwr OCMUKDWE WHY OW MNIADUPWMHOVUeENARUE Whe 


Ne 


sideevereSample printout for Ighzl4m. our ered 
frequency = 9600. 0000 mhz 
horizontal polarization 
maximum mode attenuation rate = 2.1000 db/km 
seawater temperature = 15.0000 degrees celsius 
seawater salinity = 35.0000 grams salt/kg seawater 
dielectric constant of seawater = 80. 8869 
conductivity of seawater = 4.6400 si/m 
rms surface bump height = . 2500 meters 
iflagb = 1 
air temperature = .0000 degrees celsius 
relative himidity = .0000 percent 
liquid water concentration in air = .0000 grams/meter***3 
m( .0000 ) = .0000 m-units 
gab( .0000 ) = - 0000 db/km 
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tropospheric modified refractivity profile 


m(z) 
(m-units) 


dm/dz 
(m-units/meter) 


gab(z) 
(db/km) 


d(gab)/dz 
((db/km)/meter) 


z 
(meters) 


. 0000 
. 0400 
. 1000 
. 2000 
. 3980 
. 7940 
- 5850 
. 1629 
6. 3100 
. 5890 
. 0000 
5.1190 
9.8110 
. 1190 
. 0960 


. 0000 
-3.9000 
-5. 3400 
-6. 4000 
-7. 4609 
-8. 4900 
-9.4700 

-10. 3500 
-11. 0400 
-11. 3200 
-11. 3300 
-10. 8700 
-9.7500 
-8. 8200 
-7.5600 
. 4330 -5. 8800 
), 0000 -3.6700 
. 8930 -. 8000 
- 48990 2.9200 
- 5260 7.6300 
- 5260 8. 8700 


. 00000D+00 
. 00000D+00 
. 00000D+00 
. OOOOOD+00 
. 00000D+00 -2. 
- OO000D+00 -l. 
. 00O000D+00 -. 
. OOOO00D+00 -. 
. 00000D+00 =, 
- 0000CD+00 = 
. 00000D+00 
. O0OO00D+00 
. 0O0000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
. 00000D+00 
- 00000D+00 
. 00000D+00 
. 00000D+00 


-97. 
~24, 
-10. 

25. 


5000 
0000 
6000 
3535 
6010 
2389 
5580 
2192 
0446 
0071 
~0414 
- 0762 
. 0902 
.0971 
- 1028 
~ 1075 
- 1108 
~ 11741 
- 1162 
. 1180 


. 00O00D+00 
. 0O000D+00 
. 0O0000D+00 
. OOOOOD+00 
. 00000D+00 
. OOOOOD+00 
. 00000D+00 
. 00000D+00 
. 0OO000D+090 
. O0OO00DD+00 
. OO0OCD+00 
. G0000D+00 
. 0O00COL+00 
. O0N00D+00 
. OOOO0OD+00 
. 00000D+00 
. 00G00D+00 
. OOOOOD+00 
. 00000D+00 
. OOOOOD+00 


Tetsdedede Mevevededs 


start search for modal eigenvalues 


tmesh= 5. 84278D-04 tol= 5. 84278D-05 


zeros found in expanded search rectangle defined by 
tleft= 1.87D-01 tright= 2. 80D-01 
ttop = 2.45D-02 tbot = -5. 84D-05 


qeigen( 
qeigen( 
qeigen( 
qeigen( 
qeigen( 
geigen( 
qeigen( 
geigen( 
geigen( 
qeigen( 
qeigen( 
qeigen( 
qeigen( 


wenn tt eunuu wet 


NNNNNNNHNNNHNHNE 


.89173D-01 
- 10626D-01 
- 30399D-01 
- 34888D-01 


78790D-01 
72642D-01 


. 66269D-01 


60054D-01 


-54599D-01 
-49956D-01 
~45515D-01 
- 40284D-01 
- 31896D-01 


fe oe ok 0 Oe ee ee om oo 


-02577D-02 
. 92643D-03 
. 18555D-03 
. 88232D-03 
~17520D-02 
.09210D-02 
.04444D-02 
-04991D-02 
-01015D-02 
-01089D-02 
. 69224D-03 
.71949D-03 
- 46480D-03 





qeigen( 
geigen( 
geigen( 
qeigen( 
qeizen( 
qeigent 


14) 
15) 
16) 
17) 
18) 
19) 


hou td Wout tt 


zeros found in expanded 


tleft= 9. 25D-02 

ttop = 2.45D-02 
geigen( 1) 
geigen( 2 
qeigen( 3) 
geigern( 4 
qeigen( 5) 
qeigen( 6) 
geigen( 7) 
geisen( & 
qeigen( 9) 
yeigen( 10) 
qeigen( 11) 
qeigen( 12) 
geigen( 13) 
geigen( 14) 


1) | OS | | 


zeros found in expanded 


tleft= 0. 00D+00 

ttop = 2.45D-02 
qeigen( 1) 
geigen( 2 
qeigen( 3) 
qeigen( 4) 
qeigen( 5) 
qeigen( 6) 
qeigenf 7) 
qeigen( 8) 
qeigen( 9) 
qeigen( 10) 
qeigen( 11) 
qeigen( 12) 


zeros found in expanded 


tleft= -9. 35D-02 

ttop = 2.45D-62 
geigen( 1) 
geigen( 2) 
qeigen( 3) 


tm POPIN fo ht 


. 27329D-01 
~21512D-01 


15752D-01 
C3957N-01 


.G0678N-01 
.95150D-01 


te b+ bh OF OO OO 


. 32107N-03 
»01502)-C3 
~ 79822D-03 
.643517D-02 
~9413°D-02 
.00755D-G2 


search rectangle defined by 


tright= 1. 87D-01 

tbot =  -5.84D-05 
9.94029D-02 1. 21015D-02 
1.13825D-01 1.19544D-02 
1.70951D-01 1.06612D-02 
1.83648D-01  1.04976D-02 
1.77265D-01  1.02367D-02 
1.65278D-01 1.15316D-02 
1.59615D-01 1.13678D-C2 
1.53461D-01 1. 11654D-02 
1.46761N-01 = 1.12996D-02 
1.40398D-01  1.11657D-02 
1.33415D-01  1.13914D-02 
1.27009D-01 1.22238D-02 
1.20774D-01 1.220951-02 
1.06849D-01 = 1. 20849D-02 


search rectangle defined by 


tright= 9.35D-02 

tbot = -5. 84D-05 
3. 55132D-03 1.51885D-02 
3. 84280D-02 1. 36589D-02 
9, 22631D-02 1.27898D-02 
8. 55517D-02 1. 28786D-02 
7. 80306D-02 1. 28557D-02 
7.02978D-02 1. 28209D-02 
6. 22701D-02 1. 33628D-02 
5.48544D-02 1. 36509D-02 
4. 69437D-02 1. 36532D-02 
3.01131D-02 1.42761D-02 
2. 18271D-02 1. 45585D-02 
1. 27540D-02 1. 46056D-02 


search rectangle defined by 


tright= 


tbot 


-8. 
“4, 
-5. 


0. OOD+00 


= -5. 84D-05 


35376D-02 
65663D-02 
54196D-03 


Ls 
1. 
1. 


99715D-02 
72812D-02 
54024D-02 


qeigen( 4) = 
qeigen( 5) = 
qeigen( 6) = 
qeigen( 7) = 
qeigen( 8) = 


zeros found in expanded 


tleft= -1.87D-01 

ttop = 2.45D-02 
qeigen( 1) = 
qeigen( 2) = 


zeros found in expanded 
tleft= -2. 60D-01 
ttop = 2.45D-02 


zeros found in expanded 
tleft= -3.74D-01 
ttop = 2.45D-02 


zeros found in expanded 


tleft= 2. 80D-01 

ttop = 2.45D-02 
qeigen( 1) = 
qeigen( 2) = 
qeigen( 3) = 
qeigen( 4) = 
qeigen( 5) = 
qeigen( 6) = 
qeigen( 7) = 
qeigen( 8) = 
qeigen( 9) = 
qeigen( 10) = 
qeigen( 11) = 
qeigen( 12) = 
qeigen( 13) = 


zeros found in expanded 


tleft= 3. 74D-01 
ttop = 2.45D-02 
qeigen( 1) = 


-1.54690D-02 1.57166D-02 
-2.52774D-02 1. €2835D-02 
-3.56533D-02 1. 66336D-02 
-5. 61425D-02 1. 60677D-02 
~7.03978D-02 1. §8961D-02 


search rectangle defined by 
tright= -9. 35D-02 
tbot -5.84D-05 


-9.79148D-02 
-1.14139D-01 


2. 14675D-02 
2. 38198D-02 


search rectangle defined by 
tright= -1.87D-01 
tbot -5. 84D-05 


search rectangle defined by 
tright= -2.80D-01 
thot = -5. 84D-05 


search rectangle defined by 


tright= 3. 74D-01 

thot = -5.84D-05 
3. 72978D-01 1.55999D-02 
3, 64365D-01 1.53894D-02 
3.56381D-01 1.55135D-02 
3.48768D-01 °1.52448D-02 
3.41094D-01 1.49015D-02 
3. 33811D-01 1, 46419D-02 
3. 26458D-01 1. 38474D-02 
3. 18782D-01 1. 33368D-02 
3.11377D-01 1. 31438D-02 
3.04288D-01 1. 31966D-02 
2.98015D-01 1, 32973D-02 
2.91375D-01 1. 24615D-02 
2.84946D-01 1. 22797D-02 


search rectangle defined by 
tright= 4.67D-01 
thot -5. 84D-05 


4. 61863D-01 1.91926D-02 
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qeigen( 
geigen( 
qeigen( 
geigen( 
qeigen( 
qeigen( 
qeigen( 
qeigen( 1 
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zeros found in expanded 


tleft= 


ttop = 


4.67D-01 
2.45D-02 


qeigent 
qeigen( 
geigen( 
qeigen( 
qeigen( 
qeigent 
qeigen( 
qeigen( 
qeigen( 


zeros found in expanded 


tleft= 
ttop = 


5.61D-01 
2. 45D-02 


qeigen( 
gqeigen( 
qeigen( 
qeigen( 
qeigen( 
qeigen( 
qeigen( 
qeigen( 


Huntin t waa 
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-52529D-01 
.43515D-61 
.34183D-01 
. 24640D-01 
. 15495D-01 
.06728D-01 
. 98032D-01 
. 89689D-01 
. 81495D-01 
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. 90979fP)-02 
. 86615D-02 
. 78904D-02 
. 75055D-02 
~75135D-02 
. 72964D-02 
.71405D-02 
.68946D-02 
.61861D-02 


search rectangle defined by 


tright= 5.61D-01 

tbot = -5.84D-05 
5.54181D-01 2.17228D-02 
5.43260D-01 2.13685D-02 
5.32482D-01 2.11160D-02 
5.21895D-01 2.10218D-02 
5.11€25D-01 2. 10268D-02 
5.02023D-01 2. 04408D-02 
4.91730D-01 1.97277D-02 
4.81502D-01 1, 95239D-02 
4.71618D-01 1.92732D-02 


search rectangle defined by 


tright= 6. 54D-01 

thot =  -5.84D-05 
6.44509D-01 2.47339D-02 
6.33041D-01  2.42676D-02 
6.21371D-01  2.37056D-02 
6.09727D-01 2.32642D-02 
5.98138D-01  2.29288D-02 
5.86653D-01 2. 28634D-02 
5.75869D-01  2.29331D-02 
5,65183D-01 2.22738D-02 


zeros found in expanded search rectangle defined by 


tleft= 
ttop = 


6.54D-01 
2.45D-02 


tright= 


tbot 


2 modes found on same phase line 


qeigen( 
qeigen( 


1) 
2) 


zeros found in expanded 


tleft= 


7.48D-01 


6. 
6. 


7.48D-01 


= ~5.84D-05 


68467D-01 
56153D-01 


2. 
2. 


47500D-02 
47347D-02 


search rectangle defined by 


tright= 


135 


8.41D-01 
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ttop = 


2.45D-02 


tbot = 


2 modes found on same phase line 


modal eigenvalues 


eigenvalue 


-1. 
~9. 
-8. 
nD. 
“4, 
S35 
a2. 
-1. 
-5. 


NM DN DNDN HD Be BH BR eR BR RR RR 100 1 OSI YI OD OT EO © fH WD 


14139D-01 
79148D-02 
35376D-02 


.03978D-02 


81423D-02 
65663D-02 
58533D-02 
52774D-02 
54690D-02 
54196D-03 


-55132D-03 
. 27540D-02 
. 18271D-02 
.01131D-02 
. 84280D-02 
. 69437D-02 
-48544D-02 
. 22701D-02 
. 02976D-02 
. 80306D-02 
-55517D-02 
. 22631D-02 
» 94029D-02 
. 06849D-01 
- 13825D-01 
. 20774D-01 
-27009D-01 
- 33415D-01 
-40398D-01 
-46761D-01 
-53401D-01 
-59615D-01 
-65278D-01 
. 70951D-01 
. 77265D-01 
- 83648D-01 
-89173D-01 
-95150D-01 
. 00678D-01 
-05957D-01 
. 10626D-01 
-15752D-01 
»21312D-01 
. 27329D-01 
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nrmode= 94 
theta 
.38198D-02 3.47000D-04 
.14675D-02 3.37466D-04 
.99715D-02 3.39522D-04 
.88961D-02 3.49316D-04 
.80677D-02 3.66463D-04 
.72812D-02 3.89814D-04 
.66336D-02 4.23955D-04 
.62835D-02 4. 84349D-04 
.57166D-02 5.67769D-04 
.54024D-02 7.28131D-04 
.51885D-02 9,.68348D-04 
.46056D-02 1.25460D-03 
»-45585D-02 1.53413D-03 
.42761D-02 1.76251D-03 
.36589D-02 1.96946D-03 
.36532D-02 2.16626D-03 
.36509D-02 2.33540D-03 
.33628D-02 2.48351D-03 
.28209D-02 2.63465D-03 
.28557D-02 2.77370D-03 
.28786D-02 2.90271D-03 
.27898D-02 3.01313D-03 
.21015D-02 3.12585D-03 
.20849D-02 3.23999D-03 
-19544D-02 3.34337D-03 
.22095D-02 3.44356D-03 
.22238D-02 3.53092D-03 
-13914D-02 3.61798D-03 
.11657D-02 3.71101D-03 
.12996D-02 3.79399D-03 
-11654D-02 3.87856D-03 
.13678D-02 3.95622D-03 
.15316D-02 4.02569D-03 
-06612D-02 4.09370D-03 
.02367D-02 4. 16832D-03 
-04976D-02 4.24268D-03 
.02577D-02 4.30584D-03 
.00756D-02 4. 37319D-03 
.04132D-02 4.43471D-03 
-04317D-02 4.49259D-03 
.92643D-03 4.54304D-03 
.79822D-03 4.59766D-03 
.01502D-03 4.65652D-03 
.32107D-03 4. 71922D-03 
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-5. 84D-05 


sloclestectents 
BEGETS 


. 36133D-03 
.11498D-03 


. 88036D-03 
- 64885D-03 


.41421D-03 
. 17081D-03 


. 92120D-03 
- 64625D-03 


. 35548D-03 


. 03582D-03 


. 68051D-04 


. 70059D-04 
. 64688D-04 


. 96631D-04 


. 39607D-04 
. 08626D-04 


. 86226D-04 
63474D-04 


. 38289D-04 


. 26958D-04 
17257D-04 
07852D-04 
. 89575D-04 
- 82645D-04 
75087D-04 
. 73620D-04 
.69523D-04 
.54177D-04 
.47335D-04 
.45841D-04 
- 40966D-04 
- 40704D-04 
.40268D-04 
. 27526D-04 
. 20257D-04 
. 21161D-04 
. 16655D-04 
. 12820D-04 
- 14982D-04 
. 13703D-04 
. 06994D-04 
. 37065D-05 
- 48020D-05 
. 63419D-05 





atnu db/km 


2. 
. §371D+00 
. 7091D+00 
. 6170D+00 
. 5461D+00 
.4788D+C0 
.4234D+00 
. 3935D+00 
. 3450D+00 
. 3181D+00 
. 2998D+00 
- 2499D+00 
. 2459D+00 
. 2217D+09 
. 1689D+00 
. 16864D+00 
. 1682D+00 
. 1435D+00 
. 0972D+00 
. 1001D+00 
. 1021D+00 
. 0945D+00 
. 0356D+00 
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0384D+00 


0342D+00 


. 0230D+00 
. 0448D+00 
. 0461D+00 
. 7483D-01 
~5552D-01 
.6698D-01 


5550D-01 
7282D-01 


. 8683D-01 


1235D-01 
7602D-01 


- 9835D-01 


7782D-01 


.6224D-01 
.9113D-01 
.9271D-01 
-4947D-01 
~5292D-01 
.7148D-01 
.1209D-01 
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. 30399D-01 
.31896D-01 
. 34866D-01 
.40284D-01 
-45515D-01 
.49956D-01 
-54599D-01 


60054D-01 
66269D-01 
72642D-01 


. 78790D-01 
. 84948D-01 
-91375D-01 
- 98015D-01 
. 04288D-01 
-11377D-01 
- 18782D-01 
. 26458D-01 
-33811D-01 
-41094D-01 
.48768D-01 
- 56381D-01 
- 64365D-01 
. 72978D-01 
-81493D-01 
. 8988°D-01 
- 98032D-01 
. 06728D-01 
. 15495D-01 
. 24640D-01 
. 34183D-01 
-43515D-01 


52529D-01 


- 61863D-01 
- 71618D-01 
-81502D-01 
-91730D-01 
- 02023D-01 
»11825D-01 
. 21895D-01 
. 32482D-01 
-43260D-01 
.54181D-01 
-65183D-01 


75869D-01 


- 86653D-01 
. 98138D-01 
.09727D-01 
. 21371D-01 
. 33041D-01 


- 18555D-03 
- 46480D-03 
. 85232D-03 
-71949D-03 
- 69224D-03 
.01089D-02 
.01015D-02 
. 04991D-02 
. 04444D-02 
. 09210D-02 
.17520D-02 
.22797D-02 
. 24615D-02 
. 32973D-02 
. 31966D-02 
. 31438D-02 
. 33368D-02 
. 38474D-02 
- 46419D-02 
-49015D-02 
- 524481-02 
-55135D-02 
- 53894D-02 
- 55999D-02 
- 61861D-02 
. 68946D-02 
. 71405D-02 
. 72964D-02 
~75135D-02 
- 75055D-02 
. 78904D-02 
. 86615D-02 
. 90979D-02 
-91926D-02 
- 92732D-02 
-95239D-02 
-97277D-02 
. 04408D-02 
- 10268D-02 
- 10218D-02 
. 11160D-02 
- 13685D-02 
.17228D-02 
- 22738D-02 
. 29331D-02 
. 28634D-02 
. 29288D-02 
- 32642D-02 
- 37056D-02 
- 42676D-02 
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. 75019D-03 
. 76638D-03 
. 79691D-03 
. 85181D-03 


90448D-03 


. 94869D-03 
. 99440D-03 
.04766D-03 


10756D-03 


. 16837D-03 
. 22643D-03 
. 28389D-03 
- 34313D-03 
. 40378D-03 
.46027D-03 
- 52345D-03 
- 58872D-03 
- 65564D-03 
.71907D-03 
. 78111D-03 
. 84578D-03 
- 90923D-03 


97497D-03 


. 04515D-03 
- 11381D-03 
. 18078D-03 
- 24497D-03 
. 31279)D-03 
- 38044D-03 
-45021D-03 
. 52228D-03 
. 59206D-03 
- 65873D-03 
. 72700D-03 
. 79763D-03 
- 86847D-03 
-94101D-03 
. 01331D-03 
. 08148D-03 
- 15075D-03 
- 22286D-03 
- 29558D-03 
- 36854D-03 
-44134D-03 
-51138D-03 
- 58133D-03 
- 65513D-03 
. 72892D-03 
. 80238D-03 
. 87532D-03 


- 22215D-05 
-69643D-05 
. 04647D-05 
. 80035D-05 
-67707D-05 
. 00029D-04 
. 90411D-05 
.01853D-04 
.00134D-04 
. 03472D-04 
. 10108D-04 
- 13802D-04 
. 14206D-04 
. 20498D-04 
- 18348D-04 
. 16527D-04 
. 16857D-04 
. 19895D-04 
. 25368D-04 
. 26222D-04 
. 27702D-04 
. 28557D-04 
. 26126D-04 
. 26366D-04 
- 29642D-04 
. 33851D-04 
. 34403D-04 
. 34169D-04 


34413D-04 


. 32898D-04 
. 34320D-04 
. 38626D-04 
.40447D-04 
- 39711D-04 
» 38840D-04 
-39196D-04 
. 39179D-04 
-42723D-04 
-45402D-04 
-43959D-04 
-43161D-04 
~43428D-04 
-44363D-04 
. 46577D-04 
- 49509D-04 
-47679D-04 
. 46674D-04 
-47398D-04 
-48781D-04 
. 50898D-04 
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frequency = 


9600. 0000 mhz 
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.0146D-01 
- 2439D-01 
. 7454D-01 
-46190-01 
-2945D-01 
-6509D-01 
.6445D-01 
. 9848D-01 
. 9380D-01 
. 3456D-01 
. 0057D+00 
. 0509D+00 
- 0664D+00 
. 1379D+00 
- 1293D+00 
. 1248D+00 
- 1413D+00 
- 1850D+00 
. 2530D+00 
- 2752D+00 
. 3046D+00 
. 3276D+00 
. 3170D+00 
- 3350D+00 
. 3852D+00 
-4458D+00 
- 4668D+00 
- 4802D+00 
-4988D+00 
-4981D+00 
- 5310D+00 
- 5970D+00 
. 6344D+00 
. 6425D+00 
. 6494D+00 
- 6708D+00 
. 6683D+00 
. 7493D+00 
. 7994D+00 
. 7990D+00 
. 8071D+00 
. 8287D+00 
- 8590D+00 
- 9062D+00 
- 9626D+00 
- 9566D+00 
- 9622D+00 
- 9909D+00 
. 0287D+00 
. 0768D+00 
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